KARMAN 


wry 
si. Ve 


I, S. SOKOLNIKOF! 


N. BRILLOUIN 
FELLER , 
LE CORBETILLER 
R. SEARS 


LISUNOFI 
TIMOSHENBKO 


Abii 


rAYLOR 





JULY + 1944 NuMBER 2 























QUARTERLY OF APPLIED MATHEMATICS 





Vol. II JULY, 1944 No. 2 





ON COMBINED FLEXURE AND TORSION, AND THE 
FLEXURAL BUCKLING OF A TWISTED BAR* 


BY 
J. N. GOODIER 


Cornell University 


1. Introduction. When a straight uniform slender bar is twisted, the straight 
form becomes unstable at a certain value of the twisting couple, and the center line 
of the bar becomes a space curve. Elements of the bar are bent about both principal 
axes of section, and the buckled form thus possesses strain energy of flexure as well 
as of torsion. If the bar is twisted to the critical configuration, and its end sections 
then held against further rotation, the jump to the buckled form means the appear- 
ance of flexural energy at the expense of the torsional energy. The occurrence of the 
flexure must therefore produce some relief of the torsion, that is, it must modify the 
amount of twist. 

It proves to be impossible to account for the transference of strain energy from 
that of torsion to that of flexure if the strain energy is represented in the accepted 
form of the theory of small bending and torsion of thin bars— 


1 l 
—f (Eliu’”? + EIqv'"* + GCB’*)dz, 
- 0 


where EJ,, EJs, GC are the flexural and torsional rigidities, u, v the components of 
deflection parallel to the principal axes of the section, and 6 the torsional rotation, 
as functions of the axial co-ordinate z. Coincidence of shear center and centroid is 
assumed, and secondary effects of non-uniform torsion! are disregarded, for sim- 
plicity. If for instance this form is used in the potential energy, and the differential 
equations of the bar buckled from a state of simple torsion by couples M; are found 
by means of the theorem of stationary potential energy, the correct equations? 


EI ,u"’ + My’ = 0, EIyv"’ — Myu' = 0, M; = GCp’ 


are not obtained. The terms Mv’, M3u’ in the first two fail to appear. These equations 
are nevertheless easily derived directly as conditions of equilibrium. 

The comparison with the corresponding problem of the bar under thrust is useful. 
The bar is compressed to the critical state, and the ends held against further ap- 
proach. The bar jumps over to the bent form, and energy of bending appears. But 


* Received March 24, 1944. 
1]. N. Goodier, (i) The buckling of compressed bars by torsion and flexure, Cornell University Engi- 
neering Experiment Station, Bulletin 27, 1942; (ii) Flexural torsional buckling of bars of open section, 


Bulletin 28, 1942. 
2S. Timoskenko, Theory of elastic stability, McGraw-Hill, 1936, p. 168, or (1) (ii) equations 2, 3, 7. 
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the transition to the bent form involves a lengthening of the bar, and some of the 
compressional strain energy is thus released to supply the energy of flexure. The 
Euler problem has been analysed from this point of view by R. V. Southwell.* 

This lengthening of the bar is of the second order in the derivative of the bending 
displacement with respect to the axial coordinate. It can be disregarded in writing 
down the differential equation of equilibrium, but not in energy methods. It is 
natural to look for something analogous in the torsional problem by investigating 
the nature of combined torsion and flexure to a higher order of small quantities than 
formerly. This is done in what follows and the required new terms in the strain 
energy are found. At the same time the nature of combined torsion and flexure is 
‘clarified, and the energy method is made available for more difficult problems of 
buckling from a twisted state such as those of non-uniform bars. 

2. Finite bending and torsion of a thin bar. Let the axis (of centroids) of the 
undeformed straight bar lie along the z-axis of fixed cartesian axes u, v, 2. The bar is 
now subjected to small bending and twisting. Its axis becomes a space curve, con- 
sisting of points of co-ordinates u, v, z. Even if the deflection (u, v) is small, the 
geometrical torsion of this curve is not small. The bending may be in one plane (the 
osculating plane) at one point, and in a perpendicular plane at another. 

The geometrical torsion 7, of the curve is distinct from the torsion 7 of the bar. 
When the deflection (u, v) is prescribed the space curve of centroids is definite, with 
definite curvature and torsion. The cross sec- 
tions of the bar must be in the normal planes of 
this curve, but the torsion of the bar remains 
indefinite until the orientations of the principal 
axes in these planes are specified. 

In Fig. 1 the tangent, normal and binormal 
at P are indicated by #, n, b.4 As the origin of the 
triad moves along the curve with unit speed, it 
has a component 7, of angular velocity about #, 
and a component «x (the curvature) about 3, 
right handed rotations looking along the posi- 
tive axes being reckoned positive. Define an angle y such that r.=dy/ds (s being arc 
length increasing in the sense of ¢) and y=0 at some chosen reference section s=So, 
as for instance one end of the bar. 

Let f be the angle which one principal axis p (Fig. 1) of the cross section at P 
makes with the principal normal n, positive when this axis is obtained from n by 
positive rotation about ¢. Let fo be its value at s =so. Then the rate of rotation of the 
tpq-triad about ¢ is given by t.+df/ds or (dy/ds)+(df/ds) and this is by definition 
the torsion of the bar. 

Accordingly if the bar is bent but not twisted, y+/ is a constant along the bar 
and in fact y+f=fo, or f=fo—vy. From this state we may derive a twisted form of the 





Fic, 1. 


3 Introduction to the theory of elasticity, 2nd ed., Oxford University Press, 1941, p. 443. 

‘The notation and conventions are those of C. E. Weatherburn, Differential geometry, vol. 1, 
Cambridge University Press, Cambridge 1939, p. 15. 

5 The discussion thus far, except for the introduction of the angle y, corresponds with that of A. E. H. 
Love, Mathematical theory of elasticity, 4th ed., Cambridge University Press, 1934, Ch. XVIII. The 


further development is different. 
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bent bar by introducing an angle of twist ¢ with ¢=0 at s=so, so that f=fyo—y+¢. 
The torsion of the bar is now d¢/ds. The bent and twisted form of the bar is com- 
pletely specified by the curve of centroids, which defines y, and the angle fo+@ which 
can be assigned independently. 

In the elementary theory of bending, the curvature is related to the bending 
moments by means of components along the principal axes of cross sections. If 1, ke 
denote these components along p and q (Fig. 1), we have (x can be regarded as an 


angular velocity about b) 


ki = «sin f, ke = x cosf (1) 
or 
ki = x sin (fo— ¥+ 9), ke = xk cos(fo— 7+ 9). (2) 
We have also 
7 = do/ds. (3) 
But 
k= (ull? + 9’? 4 3/2) 1/2, (4) 


primes denoting differentiation with respect to s. Also y is defined through dy/ds =r, 


and we have 


j nu’ y’ 2! 
sill ” - 
tT =ati a” of gg” i. (5) 
wl” x” gt 


With these formulas the bending and torsion of the bar are completely specified by 
the deflection (u, v as given functions of z) and the angles fy and ¢. The orientations 
of the principal normal and binormal are defined by the deflection curve, and the 
orientations of the principal axes relative to these are defined by fo and @. The 
formulas (2) and (3) may be used to specify not only the deformed state of the bar, 
but also an initial “bent and twisted” but unstressed state. The differences between 
the values of xi, ke and 7 then represent the changes of curvature and torsion to 
which the components of bending moment, and the twisting moment, will be respec- 
tively proportional. 

To illustrate this, and also the significance of fo, let the bar be circular and in a 
horizontal plane, with the principal axis p of all cross sections also in the horizontal 
plane. Then we may take for the initial state y =fp=@ =7T =k: =0, ke=x=1/r where 
r is the radius of the circle. Let each cross section now be rotated by the same angle 
a about ¢t. For the deformed state fp =a and 


kK, = r—' sina, Ke = r' cosa, r = 0. 
The changes of the components of curvature are 
r~! sin a, r—'(cos a — 1). 


When a is small, the second, the change in ke, is negligible. The bending moment 
induced is proportional to r~'a, and corresponds to m, that is, its axis is m, in the 
plane of the ring.® 

6 This problem is analysed from first principles in Timoshenko, Strength of materials, Part II, 2nd ed., 
Van Nostrand, 1941, p. 177. 
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3. Small bending and torsion of a straight bar. The formulas (2) and (3) must 
yield such expressions as d*u/dz?, d*v/dz*, dp/dz as their principal parts for small 
deformation. The object of the present investigation is to obtain terms of higher 
order as well. 

Let u’, v’, @ be small compared with 1, and let / be a suitable length such as the 
length of the bar, or the wavelength of a periodic deflection. The formulas (4) and (5) 
involve u’’, v’’, u’’’, v’’’. If the greatest absolute value of u’’’ and v’”’ is denoted 
by 7//?, u’’ and v’’ do not exceed // and u’, v’ do not exceed n, which is small. Let 
e denote the largest absolute value of 7 and ¢. Quantities not exceeding e, €//, etc., or 
quantities differing from them only by terms involving higher powers of e, will be 
denoted by O(e), O(e//), etc. 

The relation u’?+v’2+2’2=1 yields 2’? =1—O(e*?) and so s’=1—O(e’). It yields 
also 
sf = — (Wu! + 0'0")(1 — uv? — 0)? = O(e/I) (6) 


2 


and 
2’ = O(e?/I?). (7) 
Then (4) yields x= [u’’?+0’’?+0(e*//*) ]"/*. Since u’’?, v’’? are O(e2/l?) we have as an 
approximation 
k= (ull? + y!2)1/2 (8) 
in which the error is of order e?, relative to the part retained. 
The determinant of (5) yields u’’v’’’—u’’’v’’ with an error of order €*. Then (5) 


becomes 
te = (ule — lv") (ul? + o'"?)-! (9) 


with an error or order e?. 
Now the right of (9) may be identified as (d/ds) tan—! (v’’/u’’) and, in view of the 
equations defining y (dy/ds =7., y=0 when s=So) we have 


wv uw" 


9 a 
U v0 


y = tan~'—— — tan” ers (10) 
u U9 
where u¢’, v¢’ are the values of u’’, v’’ at s=so. The inverse tangents are principal 
values. The values of sin y and cos y are required. From (10) 
sot / y’)- 


tan y = (uo’v”’ — v9’ u"’) (uc uw” + 29 


and therefore 


J? 47 


4? Jf 

Vo wU 
ss ay ra 
un U 


The ambiguity of sign involved in obtaining the sine and cosine from the tangent is 
disposed of by the consideration that if v’’/u’’ slightly exceeds vg’ /ug’, both being 
positive, 7 must be a small positive angle. 

4, Expressions for small curvature and torsion. Expanding the first of (2) in the 


form 
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¢? 
kK, = c{(sin fo cos y — cos fo sin v(i _ * tae ) 


3 
+ (cos fo cos y + sin fo sin y) (+ — tee )} 


and substituting for x, cos y, sin y from (8), (11) we find 

ki = wu’ sin (fo + 4) — 0” cos (fo + 6) + o[u” cos (fo + 6) +0" sin(fo+4)] +--+ (12) 
and similarly 

ko = u’’ cos (fo + 6) + 0” sin (fo + 6) — o[u’’ sin (fo + 5) — 0’ cos (fo + 5)}+--- (13) 


where cos 6=ug! (ug’? +u9'?)-"/2, sin 6=v9! (ug’? +09 2)—!/2. In these developments 
the errors are of order e? relative to the leading terms. They are therefore accurate 
as far as explicitly carried. 

Since dz/ds = 1—O(e?), replacement of differentiation with respect to s by differ- 
entiation with respect to z, to any order, will involve errors of order e?. Thus the 
primes in the terms set out in (12) a::d (13) may be taken to indicate differentiation 
with respect to z, and the developments remain correct to this order. In the same 
way the torsion dé/ds may be replaced by d@/dz with an error of order e’. 

The angle fo, while significant of course when the axis of the bar is appreciably 
deflected, tends to become merely a rigid body rotation when the bar is nearly 
straight. In order to eliminate such a rigid-body rotation, we observe that there is 
as yet no connection between the w-axis and the principal axis p. If these axes coincide 
when the bar is undeformed, small torsion and bending, free of large rigid body rota- 
tions, will restrict the angle between them to be of the same order as ¢. Then the direc- 
tion cosines of p, relative to the u, v, z axes must be 1—O(e*), O(e), O(e) at most. 

The direction cosines of n, the principal normal, are u’’/x, v’’/x, 2’’/x so that, if n 
is the unit vector along n, i, 7, and k unit vectors along the axes of u, v and 2, 


n= «x Vy" i + vit 2k). 
The direction cosines of b, the binormal, are used as the coefficients of i, j, k in 
b = «[(v's!"” — 2'0")i + (2'u"” — u'2!"")j + (u'o” — v'u'’)k| 


where b is the unit vector along the binormal. 
Since the principal axis p (Fig. 1) is in the plane of b and n, and is derived from n 
by a rotation f towards b, the unit vector along it is given by n cos f+) sin f or 
x [u’’ cos f + (v's! — 2'v") sin fi 
+ «x [v” cos f + (2’u"’ — u's’) sin fj (14) 
+ «x [s" cos f + (u’v” — v'u’’) sin f|k 
and the coefficients of 1, 7, k give the direction cosines of p. 
The first of these is of order 1 without restriction on f. The second may be repre- 


sented as 


O(1/e) | Oce/) cos f + O(1)O(e/1) sin f — O(e)O(e?/1) sin f] : 








98 J. N. GOODIER [Vol. II, No. 2 


from which it is apparent that these direction cosines will not be small of order € unless 
’(v’’ cos f + uv” sin f) 


is small of this order. This expression may be developed, by the processes which led 

to (12) and (13) as 

x0" { u"’ cos (fo + 6) + v” sin (fo + 6) — ou” sin (fo + 6) — v” cos (fo +6) | +--+} 
+ xu!" { u!’ sin( fo + 5) — 0” cos (fo + 6) + ou” cos (fo + 6)+ v” sin (fo+6)] +--+} 


and will be small of order ¢€ only if fo+46 is small of this order. 
This result simplifies (12) and (13) to 


m= —ov'tul(otfots, c=u" +0"(o+ fot), (15) 


and with r =¢’ these constitute approximations to m4, ke, and 7 with errors of order €° 
relative to the leading terms. It is now implied of course that one principal axis (p) 
coincides with the u-axis in the undeformed state, and that in the deformation it ro- 
tates from it by an angle of the same order as uw’, v’ and @. This is the case if one sec- 
tion of the bar is fixed against rotation, or against rotation of the type ¢ only. 

The third direction cosine in (14) is of order € without further conditions. 

5. An alternative torsional co-ordinate. The angle ¢ represents a rotation of the 
cross section about #, from the torsionless configuration associated with the deflec- 
tion u, v. This torsionless state is far from being geometrically obvious, and the 
terminal values of @ and f corresponding to various types of simple end constraints 
are not immediately obtainable. 

A representation of the torsion and flexure to the second order which does not 
suffer from these disadvantages is desirable. A straight bar (initially along the z-axis, 
Fig. 2) may be imagined brought to a bent 
and twisted state by supposing it cut into 
thin discs. Let a typical disc be translated 
without rotation so that its centroid is 
brought to its final position P on the de- 





= ar 

flected curve and the principal axes are 

brought to x, 71 parallel to x, y. It must 

now be rotated so that the tangent at P to 

the deflection curve is normal to it, in accord- 

ance with the theory of flexure of thin bars. 

* Let this rotation consist of a rotation about 
ee # bringing x; to x2 in the normal plane at P, fol- 


lowed by a rotation about x2 bringing y; to ye 
in the normal plane. The configuration so produced is evidently a possible state of 
bending and torsion. The principal axis x2 is still parallel to the plane xz. This configu- 
ration is to be used as a reference from which to measure the torsional rotation of cross 
sections. To the first order the torsion is zero, but to the second order it is not. 
To determine its value, let the x, y axes in Fig. 2 correspond with the u, v axes, 
and let x2 be the principal axis p. Then, in the proposed configuration, p is everywhere 
normal to the y, or v, axis. Thus the coefficient of j in (14), which represents the direc- 
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tion cosine of » with the v axis, must vanish, so that the value of f is determined by 


the equation 
Mt 


v 
fu «= ta" 


u's!’ — 2'u!" 


(16) 





The torsion of the bar is r.-+df,/ds, rT. being given by (5), and is thus expressed in 
terms of the derivatives of u and v. When expanded in powers of these derivatives 
its leading term is u’’v’. This is an approximation to the torsion with error of order e. 
Thus if ¢; is the value of ¢ corresponding to this configuration ¢/ =u’’v’[1+0O(e) ]. 
Also, fo is obtained from (16) by putting wo, vo for u, v and it is easily found that 
tan fo= —(vi’ /ug )+O0(e). Since tan 6=v9 /ug’ it follows that fo+6=O(e*). This 
being so fy +4 in (15) ceases, for this particular configuration, to be significant, since 
its products with u’’, v’’ are of the order of the terms neglected. 

Now consider an arbitrary state of (small) flexure and torsion specified by 4, 2, @. 
It may be derived from the reference state just defined simply by rotating cross- 
sections about ¢ in order to convert ¢; to @. Let 8B be the amount of such rotation. Then 
¢—¢,=6, and r=¢’=6’+¢/, that is 

T= P+ u'r’ (17) 

with error of order e?. 

Let s now be measured from one end of the bar so that s)=0. Then @¢; like @ is 
zero at s=0 and ¢;=f-u''v'ds. Thus ¢=B+Jf-u''v'ds and the integral is of order e?. 
Moreover f)+64 is not altered by the rotation 8 so that it is still of order e?. The first 


of (15) becomes in consequence 


a 
Ko = —o' t+ u"(8 +f w'dds), 
re 
ley? 


The first term is or order €/l, u’’B is of order €?/] and u’’f\u''v'ds is of order €/l. 
Therefore, with an error of order e® the new formulas for the components of curvature 
are 

kK, = —v’ + Bu”, ko = ul’ + Bo". (18) 
These with (17) give an alternative representation of the torsion and flexure, con- 
venient because fp and 6 have been eliminated, and 8 is relatively easily envisaged- 
being the angle by which the cross section must be rotated, about the deflected, tan- 
gent, to bring p from the position parallel to the axial plane in which it originally 
lies, to its final position. At fixed ends 8 is clearly zero. 

6. Energy considerations. The strain energy is given by 


1 I 
—{ (EI yx + El ox? + GCr*)dz. (19) 
2/0 


The integration with respect to z rather than s will involve an error of order e’. 
Consider now the problem referred to in the introduction—the straight bar twisted 
until it buckles. Let the state just prior to buckling be 


68 = B, = 0, v= 0, 


and after buckling 


B= B+ 6, “= 1, v= 2}. 
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Then B is small in the sense of ¢@ in the preceding analysis. But 6), 1, v1 are to be 
true infinitesimals, since we seek a buckled form which comes to the straight form as 
a limit. Thus they are to approach zero after a fixed value has been assigned to B. 
The expressions (17) and (18) are now used in (19), and terms to the second order 
in %, 2%, 6; and their derivatives, without regard to B, are retained. The result is 


1 l 
—{ [El,(ui ’ oo 2Buj' vz") 4 EI.,(v;'? — 2Bu{' v{') 
- 0 
+ GC(B’? + 2B’B/ + Bi? + 2Buf'v{) |dz. (20) 


Let M; be the critical torsional couple GCB’. On buckling, some work is done by 
this couple, but exactly how much, in terms of §;, 1, 1: depends on the end constraints 
of the bar. 

If the ends are in bearings which constrain the axis of the bar to remain fixed in 
direction at the ends—i.e., the ends are “built-in” with respect to flexure—the rotation 
of one end may be set as zero, and that of the other about the axis is then the value of 
B, at that end. The potential energy of M; in the buckled form is — M;f.B{ dz referred 
to the twisted but unbuckled form as zero. The total potential energy is thus this 
term together with (20), omitting + ['GCB! *dz which is the energy of the unbuckled 
twisted form. 

If the potential energy is now varied by varying ™ to u,+€.m(z) the coefficient 


of €, in the variation of the potential energy is 
f [— El.Bo{’ + El,(u{’ + Bo{’) + GCB'v{ |ni dz 
and this must vanish if the buckled state is a possible state of equilibrium. Since 


Bv{’ is small compared with uj’, on account of the smallness of B, the conclusion 


is that the equation 


EI uj’ M301 = 0 (21) 
must be satisfied. Similarly variation of 1 yields 
ETI ov{' _ M3ui == (). (22) 


Variation of 8, yields GCB’+GCB/ —M;=0, that is Bf =0. Equations (21) and 
(22) are identical with the equations obtainable by direct equilibrium considerations. 
They are derived in this manner here in order to show that the terms M;v/ , — M3u/ 
arise from terms in the strain energy of torsion which are of higher order than the 
term 3/,GC8’*dz hitherto accepted. It is to be expected therefore that in (17) and (18) 
the terms of the second order will be requiréd in energy calculations in other prob- 
lems where torsional loads cause, or contribute to, buckling. 

When the equilibrium of the straight twisted form is neutral, the work done by M; 


during buckling is equal to the gain of strain energy. Then 


l 1 j 
us [ Bidz = . f [ET,(ui{'? + 2Bu{’o{') + El.(v{'? — 2But' v1‘) 
+ GC(2B'Bj + Bi? + 2B’uf'v{)|dz. (23) 
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The term Bus’ vj’ in the flexural terms is small compared with uj” or v;’ and will be 
dropped. Introducing M;=GCB’ the resulting equation yields 
“ 1 fi(Eliui!? + Eli’? + GCB{*)dz - 
as 2 Shut! of dz { 

Now equations (21), (22) (after one differentiation) together with B/ =0 are the 
Euler differential equations for the functions ™, 1, 6: making the right of (24) a 
minimum. Since B/ =0 the term GC/? in the numerator of (24) may be dropped. 
The critical M; is the least value of the right of (24) with or without this term. With- 
out it the equation may be interpreted as showing that the energy of flexure which 
appears when buckling occurs is accounted for by a decrease of torsional energy of 
amount M,f-ui' vi dz. 

The same equation is suitable for the approximate determination of the critical 
torque by the Rayleigh method—assuming simple plausible forms for u; and », and 
adjusting the parameters of these forms to obtain a least value of M3. This method 
is applicable to non-uniform bars. 

Equation (23) would in general require modification if the ends are not “built-in,” 
for instance if they are attached to Hooke’s joints. For then the work of M; is not 
done merely on a rotation Spi dz. Certain terms of higher order must be added to fy, 
and these can be of the same order as uj‘ vj. Such terms would be significant in (24). 
Nevertheless (24) is appropriate in the Rayleigh method whatever the end con- 
straints, for its minimizing conditions are the differential equations of equilibrium 
which must be satisfied irrespective of end constraints. 

There are expressions other than the right of (24) which yield the critical M; as 
a minimum value. If (21) and (22) are multiplied respectively by u;" , o{’ , integrated 
along the bar, and added, the result yields another in the form 


Si(Elyul!? + Elvi"?)dz 


3= 
Shut of! — vi uj’ )dz 
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MEMBRANE STRESSES IN SHELLS OF CONSTANT SLOPE* 


BY 
VLADIMIR MORKOVIN! 


Brown University 
. 


1. A surface S of constant slope may be generated by a straight line LZ sliding 
along a plane curve C> (say, in the xy plane), maintaining a right angle with the tan- 
gent to Cy and a constant angle @ with its binormal (i.e., with the z axis). When a 
closed curve Cp is chosen, the surface is an obvious generalization of a circular cone? 

(see Fig. 1). Since “near-conical” shells occur 

. Z often in practice,’ it may be of interest to dis- 

cuss such effects as fall within the scope of the 
membrane theory of sheils. 

We introduce the following notations: 

i, j, k, unit vectors in fixed rectangular 

directions x, y, 2; 

A, A, FB unit tangent, normal, and binor- 
mal of curve Co; 

t, length along generators L; 

Sty Pty arc length and radius of curva- 
ture of a horizontal section C; of 
the surface S; subscripts 0 and 
1 will designate corresponding 
quantities in the end sections Co 
and C, of the shell; 

#=7(s9), vector equation of curve Co; 

Y, angle between the positive x axis 
and the outward normal of Co; 

E,v, G, Young’s modulus, Poisson’s ratio, 











and shear modulus; 





h, thickness of shell having the sur- 
face S for middle surface; 
Fic. 1. N,, Ni, normal forces per unit length of 


sections of the shell which are per- 
pendicular to s- and ¢-directions respectively (Fig. 3); 
Naty shearing force in s-direction per unit length of shell section perpendicular 
to t-direction; 


* Received Oct. 16, 1943. 

1 The author wishes to express his appreciation to Professor W. Prager for proposing the problem 
and for other valuable suggestions. 

2 Non-circular cones (for which the generators meet in one point while their “slope” varies) have 
been considered recently by A. Pfliiger, Z. angew. Math. Mech. 22, 99-116 (1942). 

’ The fuselages of some aeroplanes, for instance, can be approximated by one or several shells of 
different slopes connected by stiff bulkheads. The construction of models is relatively simple because each 


portion forms a developable surface. 
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Ces, Crt, Cst, Strains corresponding to N,, N,, and N,:, respectively. 
We note some simple relationships: 


a (1.1) 
rt =}; p= R; ? 
X= —ising+jcos¢; fi = —icos¢ —jsing. (1.2) 


Since po=dso/dg, we obtain from (1.2) the Frenet-Serret formulae for a plane 
curve: 


dX/dso = Z/po; dji/dso = X/ po. (1.3) 
The vector equation of the surface of constant slope S has the form: 
R(so, 2) = #(so) + t(@ sin 6 + 9 cos 8). (1.4) 


For a constant value of t, (1.4) is the vector equation of the horizontal section C;. 
Then, 0R/ds, is the unit vector tangent to C;. Since 0R/ds,=K(po—t sin 0)dso/pods:, 
C, is parallel to Cp at corresponding points (see Fig. 2), and 


ds,/dso = (po — ¢sin 6) /po. ; (1.5) 





Fic. 3. 


Hence, for corresponding points, the centers of curvature of Co and C, coincide, and 
Pt = po — tsin @. (1.6) 


If the shell is long, it may happen that at some point p;=0. At such a point the tan- 
gent to C; ceases to turn continuously (see points P, P’ in Fig. 2). We shall discuss 
only the portion of the shell where ¢ sin @<po, i.e., the open shell without the “tail 
edge.” 

2. An element of a shell of thickness h having the surface S for middle surface is 
shown in Fig. 3. According to the usual assumptions of the membrane theory of 
shells,‘ the bending stresses as well as effects of curvature of S are disregarded and 


4 See for instance S. P. Timoshenko, Theory of plates and shells, McGraw-Hill Co., New York, 1940, 
p. 356; also the first chapter of W. Fliigge’s Statik und Dynamik der Schalen, J. Springer, Berlin, 1934. 
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one has N,,=N,,. The total forces acting on the faces fds; and hdt of the element are 


respectively: Y 
— {N.(m sin 6 + 5 cos 6) + Nad} (po — ¢ sin 0)dg, (2. 1a) 


— {N+ Nala sin 6 + 9 cos 6) } dé. (2. 1b) 


Let P=P,X+P, (g sin 6+9 cos 0)+P,(@ cos @—5 sin 0) represent the load per unit 
area of the surface. Then the condition of equilibrium of the element of the shell is: 


F) 
= { [.(a sin 6 + 9 cos 0) + Nasd](p0 — ¢ sin 0) } dtdy 


0 é = 
+ oe { N.X + N,.(a sin 0 + 3 cos @) | dtdy = (po — ¢ sin 0) Pdidg. (2.2) 
gy 


Equating the components of these forces in the n, s, and ¢ directions, we obtain three 
equations for the determination of the three stress components: 


N, = (po — ésin 0)P, sec 8, 


: { Nii in )} —N 0 = ( 6)P. ec 

— }Nai(po — ésin —N,; sin @ = — ¢sin ._- 

Ot _ ia dg (2.3) 
0 ONs 





a {N1(p0 — tsin 6) } = — + (po — ¢sin 0)P; —N, sin @. 


Ot 
We proceed to solve equations, (2.3) with the simplifying assumption that the load P 
does not vary along the generators L, and obtain:5 
N, = (po — é sin @)P, sec 0, 
f(¢) sin 0 


Nat = —————— — 3 csc (po — ¢ sin 0)(P, — P,i sec 60) + 3p¢ P, csc 8 sec 8, 
(po — # sin 6)? 


— 1 (¢) 4 
Ni — | - ste) | (2.4) 


po — tsin@ Lpo — ésin#@ 





t csc 6 = =e _ 
ang [be Ps — Sd Pr sec 6 — 4pc"P,, sec 6] 
Po — ¢81n 


— 4 csc 0(po — ¢ sin 6)[P, — P, tan @ — 4P,” csc 6 sec 0+ 4P! csc 6], 


where f(g) and g(v) are arbitrary functions of g and the prime denotes differentia- 
tion with respect to ¢. If the curve Cp is closed the continuity of stresses demands 


that f and g’ have a period of 27. 
When the load on the shell is applied only through the end sections Cy and C; the 


stress system becomes: 
sin 6 —1 ; 
N, = 0; Nu = ot > N; = 7 | f 7 _ :| . (2.5) 
(po — ¢ sin 6)? (po — ¢sin 0) Lpo — é¢sin@ 
Substituting (2.5) into (2.1a) and integrating between 0 and 27, we obtain the re- 
sultant force F; acting on the section C,; the expression for F; simplifies readily by 
virtue of (1.2): 





5 In the case of cylindrical surfaces, @=0 and integration of (2.3) leads to the special solution: 
N.=pPn; Nu =f(s)+t(P;—dN,/ds); Ni=g(s) —tdf/ds +tP,+0/2(d*N,/ds*? —dP,/ds); where f(s) and g(s) 
are arbitrary functions of the arc length s. In this connection see pp. 66-76 of Fliigge’s book. 
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27 , 
r= [ {—|—~— Gasina + » cos) |] + g(a sino + vost) bay 


po — tsin@ 
= — sin oil g’ cospde +s f_ g’ sin ede\ aa k cos 0{ g(2x) — g(0)}. (2.6) 
0 0 


The resultant moment M, about the origin due to the forces on the section C, is 


found similarly: 


M, = f R x {- |—..-, (Z sin @ + Dcos | + g'(#sin@+ rcos bdo, 
0 


po — tsin@ 


It follows by integration by parts that 


2r 2r 
M. = R00) x Ft kein f fae + cos f (icos yg + jsin g)f dg 
0 0 


_ fox x { {ve sin 6+ 3 cos ade (po — # sin @)dy. (2.7) 
0 0 


The results (2.6) and (2.7) will form the basis of analysis in later sections. 
3. Let the vector of infinitesimal displacement be 


D = ud + o(f sin 6 + ¥ cos 6) + w(f cos 6 — BP sin 8). (3.1) 


The strains in the surface are given by the following scalar products between the rates 
of change of the displacement D and the unit vectors in the ¢ and s directions: 


aR aD aR aD dR aD AR oD 
+. = —--——> Coo = ——-——? Cet eager ° (3.2) 
ot at OS; OS: Os; Ot “at “OS: 
We evaluate (3.2) and substitute the results into Hooke’s Law: 
1 Ov 
oe th ee eS 
Eh ot 
1 1 
— {N, — »N,} = ———— {w’ — (vsin 0 + wcos6)}, (3.3) 
Eh po — ¢tsin@ 
2(1 + v) 1 { Ou 
—_—— N,, = —————_- — tsiné)—+usiné+0'>. 
a" an a 


Equations (3.3) are easily integrated to yield expressions for the displacements: 





1 t 
y= fow — yN,)dt + A(g), 








Eh 
2(1 + 
“= mA “a —tsing fi eet 
oe (3.4) 
—tsiné@) rt {Ni — vN!)dt F 
| i= - dt — A'(g) csc 0 + (po — ¢ sin 6)B(¢), 
— ésin 6)? 


Il 


w= u'secO0+vtané — — (po — # sin 6)(N, — vN,) sec 8, 
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where A (v) and B(¢) are arbitrary functions. When the stresses have the form (2.5), 
the displacements can be expressed directly in terms of the functions f and g: 





csc 6 : — : 
v = — {= f'pr! + dfo'nr? — g’ Inps} + A, 

Eh 

csc? 6 : tices ; 
OT EI {(1 + v)f pr! sin? 6 + 3f"or! — 3(fp” + 3f’0’)pr? + ther? 

Ht 

+ 42'p’pr! + g"(In p: + 1)} — A’ csc 0+ pB, (3.5) 

sec 6 csc’ 0 a . 3: 
_— E} {sin? 0[— $fp’pr? + 2f’or! + g'(In pe + v)] + 3/07 

LH 


t. $(fp” + 4f’p " + 6f"p \ pr? + ps (15/’p 12 + 10fp’p”) pr? 
= 3 fo" pr 4 + 1(g'p ” + 32" ‘pr! ca Lo! p!2p7? + ge” (In pt + 1)} 
— tan 0(A + A” csc? 6) + B’p; sec 6 + Bp’ sec 0. 


Expressions for displacements D., D,, Dz in the x, y, 2 (or any other) directions are 
best derived by taking a scalar product between a unit vector in the given direction 
and D of (3.1). For instance, 

.=k-D = v0c0s 0 — wsin8. (3.6) 


4. The current literature on shells contains very little on the boundary conditions 
in the membrane theory of shells. We recall that local bending of the shell was dis- 
regarded according to the simplifying assumptions of the theory. Thus we cannot 
expect to satisfy all of the usual boundary conditions. For instance, we cannot ask 
that the heavy end bulkhead be considered rigid; in bending of the shell as a whole 
this would entail e,,=N;=0 in the end section which could consequently transmit 
no bending moment. By allowing deformations in the plane of the end sections we 
remove the restriction on N; and the problem of bending has a solution (see section 5). 
One has to decide in every particular problem which boundary conditions correspond 
more nearly to the assumption of no local bending. 

A casual reader might be tempted to interpret the contribution of A and B to 
the displacements in (3.4) as that of rigid body motion since it is present when the 
stresses vanish. However, it is conceivable that a given state of stress induces inex- 
tensional displacements other than those of a rigid body as necessitated by the shape 
of the shell. Thus, in the case of a non-circular cylindrical shell under torsion, A ac- 
counts for the warping of the cross-sections. ® 

In general, these inextensional deformations are accompanied by local bending 
stresses which must be small to be neglected in accordance with our assumptions. One 
would expect that no energy is expended in the inextensional deformations. The strain 
energy in shells loaded through the end-sections is 


9 
4 
N 


Ni ] 
do dt, 4.1 
v-— fof (gt +=) ode on 


6 Specifically, A =(T'/2A oGh) | [? pde—(1/2A 0) f(xp cos p+yp sin ¢)dy}. This expression is found by 
0 0 
the method indicated in section 8. 
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or 
t=t) 


2r 
V= ~f (Nat + Nw)pide : (4.2) 
2/0 | t=0 
Substituting (2.5) and the contribution due to A and B into (4.2) and integrating by 
parts, we find that our expectation is verified. The accompanying local bending, how- 
ever, absorbs energy and, therefore, places limitations on the inextensional displace- 
ments according to the principle of minimum strain energy. The minimum expendi- 
ture of energy in bending occurs when the inextensional displacements reduce to rigid 
body displacements. 
One can easily verify that the most general functions A and B corresponding to 
rigid body displacements have the form 


A = — a,sin @ cos gy — ay sin @ sin g + a, cos 0 + azo Cos 0 
— ayXo cos 8 + a, sin 0(yo cos gy — xo sin ¢), (4.3) 
B = a, cot @ cos ¢+ a, cot @ sin g + a,, 


where dz, ay, @: represent the infinitesimal translations in the x, y, 2 directions; 
Q@z, Q,, a, the infinitesimal rotations about the x, y, 2 axes; and xo, yo the coordinates 
in the base section Co. 

Instead of imposing conditions on the displacements, one may prescribe a sensible 
distribution of stresses at the boundary. We note that by (2.5) the state of stress in 
the whole shell is determined as soon as the stresses N; and N,,; are given at one end- 
section. Thus two different stress distributions which are statically equivalent over 
an end-section will determine distinctly different stress distributions in the rest of 
the shell.’ ; 

5. We shall study first the effects of taper* as exhibited in a conical shell of circu- 
lar cross-section; later, we shall discuss the influence of a variable radius of curvature 
p: of the section C;. 

Let M represent the bending moment (causing tension for x,>0) applied to the 
shell through the end-sections Cy and C,. We shall try to satisfy the conditions that 
the end-sections (bulkheads) remain plane, i.e., 


D, = 0 for ¢ = 0: D, = B(xo — t; sin @) for ¢ = t; (5.1) 


where £ is the (undetermined) angle of bending, and that the displacements due to A 
and B reduce to rigid body displacements (4.3). By virtue of (3.6), (3.5), and (4.3), 
we obtain for the first of conditions (5.1) 

sec 6 csc 0 f 1 2f(d + sin? ) +f") 
sw ee es ( sin” 

Eh \2 ; 
+ (¢’ + g’”) Inr + vg’ sin? 6+ gn} — ar cos g = 0. (5.2) 

’ This is the price that has to be paid for the simplifications due to the assumptions of the membrane 
theory. A “disturbance” of the state of stress on one end-section (the difference between the equivalent 
stress distributions) “propagates” itself along the generators without “dying out.” The general theory of 
thin shells would lead to differential equations of higher order; for these one can find solutions representing 
disturbances that die out with the distance from the end-section. 

8 All the results of sections 5=9 simplify to the corresponding expressions for a cylinder as the taper 


approaches zero. 
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By symmetry, the functions g’ and f’ are odd in x9; let their Fourier expansions read 
g’ = D> (2m + 1)aang1 cos (2n + 1)y, f’ = D> (2m + 1)bang1 cos (2n + Ie. (5.3) 
0 0 


Since the resultant force Fo on Co must vanish, one concludes from (2.6) that a,=0. 
Equation (2.7) yields M,=jrb; cos 6 or bs= —(1/2)M sec 0. It follows from the co- 
efficient of cos g in (5.2) that 
ay = (1/27Ehr*) sec? 6 csc 6(1 + 2 sin? 6)M. (5.4) 

Substituting (5.4) and (5.3) into the second of conditions (5.1) and equating coeffi- 
cients of cos g in the two members, we obtain 
M sin 0(2 + csc? 2 1 1 \ 

2rEh cos? 6 


(yr —t, sin 8)? r’? 
For the coefficients d2n41 and ben41, m >0, one obtains a system of two homogeneous 
equations with a non-vanishing determinant. Therefore don41=Denyi1=0, n >0, and 


= 6.5 














M sec 6 cos 9 M tan @ sin 9 
N i= r ’ d I et =e : ’ (3: 6) 
n(r — ¢ sin 6)? a(r — ¢sin 6)? 
M sin 0(2 + csc? 0 1 1 
2rEh cos? 6 (r — ¢sin 8)? y? 


If we designate by J; the moment of inertia, r(r —¢ sin 0)’, of C; about the neutral axis, 
we can write N,=(1/J,) Mx; sec 0. Essentially the stresses in the z direction fol- 
low the classical beam formula; the influence of taper is manifested by the presence 
of the x components of stresses N; which have to be balanced by N,:. From (5.7) we 
see that all sections C; remain plane. The rate of change of the angle of bending 
increases as the shell grows'narrower: 




















dp M(1i + 2 sin? 6) M(1 + 2 sin? 6) (5.8) 
—— = . J 
dz mEhcos* 6(r — ¢ sin 6)* EhI, cos* 6 
Further effects of taper are apparent in the other displacements: 
M tan @ 2 csc? 6 1 ; 
v = —— cos |< __ - — aeato— nh, (5.9) 
2rEh r — tsin 6 r 
M sec 0 csc? @ — 2 — 2p 1 
w= sin of ; — — (2 csc? 0 — ») 
2rEh r—tsin@ r ' 
+—( — tain 6) (csc? 6 + ak, (5.10) 
y2 
M sec? 6 f csc? 6 — 4 
w=- - cos g< — — — — cos? 6(2 csc? 0 — v) 
2rEh lr — tsin@ r ' 
+ —(r — tsin 6)(csc? 6 + nt, (5.11) 
r 
M sec@ (2 —csc?6@+ 2vsin?g 1 
D, = \-— + — (2 csc? 6 — v) 
2rEh r —tsin@ r 


1 
~— (r — tsin 6)(csc? 6 + 2. (5.12) 
r 
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If we take for X, the fictitious displacement of the axis of the cone, the average of D, 
over C,; (by analogy with a cylinder or prism), we obtain for the slope of the deformed 
axis 





dX M sin 6 ‘eee 


1 
— = — (2+ csc? @)?. 5.8d 
dz 2xEh cos? 6\ (r — ¢sin 6)? " r? + yt ( 


Comparison with equation (5.7) shows that the axis is not perpendicular to the sec- 
tions C; as one might expect. Nor is the increment in slope equal to 8, the angle be- 
tween the end sections. In fact, for csc? 0@=2+-, a large taper, the axis remains alto- 
gether straight despite the angle between Cp and C,. This is due to a slipping effect 
caused by an interplay of the shearing forces N,; and the x components of N,. Finally, 
let us check (5.5) by the customary® application of Castigliano’s Principle, dV/08M =B. 
Substituting (5.6), (5.9), and (5.10) into (4.2), we have 








| M?*sin@ 2p (5.14) 
aula 4rEh cos? 6 r? j 
M? sin 6(2 + csc? @ + 2») 1 1 
V= { - - -} , (5.15) 
4rEh cos?.6 (ry — 4, sin 6)? r’? 
M sin 6(2 + csc? 6 + 27) 1 1 
By = { . -—>. (5.16) 
2rEh cos? 0 (ry —%, sin 0)? r’? 


The discrepancy between (5.5) and (5.16) is negligible in practical applications, but 
is interesting theoretically. It springs from a loose interpretation of Castigliano’s 
Principle above, which is strictly true only for a concentrated couple M. Since M is 
distributed over the end sections, it does work not only in bending the shell but 
also in deforming the end-sections within their planes, as seen from (5.14). When 
the end-sections are alike as in a cylinder or prism, as much energy is spent in the def- 
ormation of one end as is gained at the other end; then, Castigliano’s Principle holds 
even for a distributed moment. But to obtain the correct angle of bending in the case 
of a cone, one must deduct from the total strain energy (5.15) the net energy absorbed 
in the-plane deformation of Cy and C;, namely 





My sin 0 { 1 1 } 
2nEh cos? @ \(r — t; sin 6)? 2 J 


6. We derive easily the expressions for the stresses in a cone twisted by a torque T 
by making either D,=0 or N,=0 at Cy and C; and using (2.6) and (2.7), 


T 
Nae = N; = 0. 6.1 
sii 2x(r — ¢ sin 6)? ; ; vail 





From the displacements or the strain energy we obtain the total angle of twist y and 
the angle of twist per unit length of the cone 


® See for instance Timoshenko, Strength of materials, vol. 1, D. Van Nostrand, New York, 1940, 
os. S92. 
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9 


ee ee ee ee “i oer apa 
r—t,sin 6)? Fr dz 2xGh(r — tsin 6)’ 


4nGh 


Here, the effects of taper as manifested in (6.1) and (6.2) are not unexpected. 

More interesting is the case of a cone supported at Cy and bent by a force R (in 
the x direction) distributed over C;. We learn from (2.6) that the function f(g) and 
hence the shear stress N,; do not actually contribute to the resultant R acting on 
any section C;. Expressions (2.6) and (2.7) show that the term in cos ¢ of g’ alone in- 
fluences the resultant force as well as moment on C;. We superpose a state of stress 
given by (5.6) with M=R cot @(r—4, sin @) in order to bring the moment across C, to 
zero, and obtain the final result 


T csc 6 1 1 dy T sec 0 
= ‘ \ ear ° (6.2) 


. — R(t; — t) cose — R(r — t; sin @) sing 
‘a enema coe aaa na oon i icin oi Sn Salli 


: —______—___— (6.3) 
(yr — t sin 0)? n(r — ¢t sin 6)? 





7. Let us now consider shells with non-circular cross-sections C;. The coordinates 
of points on C; are expressed in terms of p; and ¢ 


? ¢ 
x, = x,(0) — f pi sin ¢ dg, v= f p: cos g dg. (7.1) 
0 0 


It is clear from (7.1) that p; cannot contain any terms in cos ¢ or sin ¢ if the shell is 
closed. If only cosine terms appear in the Fourier expansion 


(7.2) 


oo 
P= %— y fT, COS Ng, 
9 


the section C; is symmetric with respect to the x axis. The simple section, for which 
r,=0 if n+3, approximates the cross-section of many a fuselage: 


r, cos g + }r3 cos 2g — 313 cos 49; 


x = 
P 3 r3 
x,(0) = 7; fe — 9 x1(71) =-—T?; oe —; R 
8 8 (7.3) 
¥: = 7, sin g — fr3 sin 29 — 31; sin 4¢; yi(w/2) = ry. ’ 


The neutral axis of the sections coincides with the y axis (i.e., is independent of ¢) 
if x, contains no constant term and if 


Pane n(Tn41 = Te 1) 


f topode = 2 8 ——___————— = 0. (7.4) 
0 


2 nN 


In bending, only sections satisfying (7.4) will be considered. 

8. The stresses in a shell of constant slope under torsion are determined from the 
conditions that the load is applied in such a manner that only shearing stresses 
are generated at the end-sections. The conditions N,=0 at t=0 and t=t, yield 
f =kpo(po—t, sin 8) and g=k(po—t, sin 8). Substituting the expression for f into N.:, we 
find the torque T on Co 


2 a aa dR or - 
T = f R X Nudpode = & sin Ag R X — dso — ; sin of RX hag} , 
0 dso 0 
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which reduces to 
T = ksin 0{2A9 — tLe sin 0} = k sin {Ao + Ai — afi sin? 6}. (8.1) 


Here the A’s represent the areas of the sections and Lp is the length of Co, all quanti- 
ties easily measurable. Then, 














T pop 
Nu =—— ae (8.2) 
(Ao + Ai — wf sin? 0) pp? 
—T t P1 4 — T sin@ t(t, _ t)p’ 
N; =— oe =(*) - ° , < (8.3) 
(Ao + Ai — 7 sin? 6) ps \ px (Ao + Ai — af? sin? 6) p> 


The effect of the variable radius of curvature of C; is observed in the expression for N;; 
tensile stresses increase directly with p’ and inversely with p?. 
The expression (4.2) for strain energy takes the form 


V were f (1+ v) sin? 0(L>) + Ly) + i sin? 8 f( a 
= —¢ é;( ») sin B — —_— 
2Eh \- yt feo -. fc har 


i, sin 0 f2* 1 1 id pi 
-.. =f (— ‘i -) ie xs f p’? In (“) io} , (8.4) 
2 0 Pi Po 0 Po 


and the angle of twist is 
T ae of "(8 sin? 6 (- *) 
y= - —_—___———| { v CSC‘ — OE eee 
Eh(Ao+A1—7#? sin? 6)? . ; 0 12 pz pa 


ty sin 0 r 1 1 jas Pl 
- p (—+—) —p’* In (=) hae. (8.5) 
2 Pi Po Po 


The quantity in the braces is of the order of sin’ @. Also, each of its terms contains the 
factor p’*. In the common case of small taper and nearly circular shell we may use as 


th 








10 








a good approximation 


Tty(Lo + Li) 


° 8.6 
2Gh(A 9 + Ay —_ nt sin? 6)? ( ) 





Yapp _ 


Neglecting the terms in the braces of (8.5) is equivalent to disregarding the effect of 
the stress NV; see (4.1). 

The inextensional displacements given by A and B in (3.5) can be determined 
from the twist of the end-sections (centers of twist along z axis) 


a=0, ¢=0; u=y(x,cosg+ ysing), t=h. (8.7) 


These displacements include warping.!® The actual process of solving (8.7) is quite 
tedious even when a definite section is given. 
9. We conclude with a short discussion of stresses in a general shell of constant 
slope bent by couples M as in section 5. We assume that the moments at the end-sec- 
10 For a treatment of warping along similar lines see R. V. Southwell, On the torsion of conical shells, 
Proc. Royal Soc. London, (A)163, 337-355 (1937). 
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tions Cy and C; are applied in such a manner that the stress N; at these sections is 
proportional to the distance from the neutral axis: 


Ni: = eo%0 for ¢ = 0; Ni = e:%, for ¢ = hk. (9.1) 
Conditions (9.1) and the fact that the moments across Cp and C, are alike lead us to 
the following expressions: 


M pop {c | —— | poo ai ae : (9.2) 


= Eudes Io I; 


where the J’s are the moments of inertia about the neutral axis of the full respective 
sections and Q, = frxpdep the variable first moment (about the same axis) of the sec- 
tion included between 0 and ¢. The expressions for the stresses themselves read: 





t; sin 6 cos 0 





Ne = ae Ee “t, (9.3) 
tip?cos8 (Ip ok 
aa M(t, — =P : 4 
t; cos 0 p? (To I; 
M(t, — é)x0 pi Mt xpi 


(9.4) 











thIgcos@ p? tt, cos@ Ip? 
t t 


The corresponding expressions for strain energy and displacements are very cumber- 
some and can hardly be useful in practical applications. 
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NON-HOMOGENEOUS STRESSES IN VISCO-ELASTIC MEDIA* 


BY 


T. ALFREY 
Monsanto Chemical Company 


1. Introduction. The purpose of this paper is the extension of the theory of elas- 
ticity to include visco-elastic media. The materials considered in this paper are 
istropic and incompressible, and are characterized by linear relations between the 
components of stress, strain, and their derivatives with respect to time. As in the 
classical theory of elasticity, only small strains will be considered. Body forces, in 
particular inertia forces, will be neglected. 

In the following, ou (t,k =1, 2,3) and €% denote the components of the tensors of 
stress and strain with respect to a system of rectangular axes xj. Ou, O22, 033 are the 
normal stresses, ¢1=021, 023=032, 031=013 the shearing stresses. Similarly, €1, €22, €33 
are the normal strains, €12 = €21, €23 = €32, €31 = €13 the shearing strains. If u; are the com- 
ponents of the displacement vector, 


é, = 3 (i,k + Uk,i)s (1) 


where the index after a comma denotes differentiation with respect to the correspond- 
ing coordinate x, i.e., u:,4=Ou:/Oxx; Uk, =OUuz/OX;. 
Irrespective of the mechanical properties of the material, the stresses must satisfy 


the equilibrium conditions 
Cikk = 0, (2) 


where the summation convention of tensor calculus has been used. Similarly, the 
strain components must satisfy the conditions of compatibility, 


€ik,lm + Elm,ik = €il,km + €km,ils (3) 


where €ix,1m =07€;,/0x,0xm, etc. While there are obviously three equations of equilib- 
rium (corresponding to the three values which the subscript ¢ in (2) can assume), it 
may at first glance appear that there are 34 equations of compatibility. On account 
of the high degree of symmetry in (3), the number of equations of compatibility re- 
duces, however, to six; three equations of the type obtained from (3) when e.g., 
i=k=1 and /=m=2, and three equations of the type obtained from (3) when e.g., 
¢=k=1, /=2, m=3. 

By themselves, Eqs. (2) and (3) are not sufficient to determine the states of stress 
and strain in a body subject to given surface stresses. A further necessary set of equa- 
tions are those relating the stress components to the strain components in the general 


* Received Dec. 21, 1943. This paper was presented at the meeting of the Society of Rheology on 
October 29, 1943. The author wishes to express his gratitude to Dr. W. Prager for his help in the prepara- 
tion of the present manuscript. Much of the material presented in this paper will appear, in modified form, 
in the book: The mechanical behavior of high polymers. By T. Alfrey. Interscience Publishing Co. 1944. 

1 According to this convention o4%,, stands for the sum of all the terms obtained by giving k the 
values 1, 2, 3. In general, whenever a subscript appears twice in the same monomial, this subscript is to 
be given the values 1, 2, 3 and the resulting terms are to be added. Such a repeated subscript is called a 
dummy subscript. 
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case of combined stresses. It is through these stress-strain relations that the properties 
of the material enter the problem. 

In the case of an incompressible material, €;; = €:;-+ €22-+€33 =0. The stress-strain 
relations are most easily discussed when the following decomposition of the stress 
tensor is introduced. Define the mean normal stress as 


o = 30% = 3(o11 + o22 + 33), (4) 
and the deviatoric part of the stress tensor as . 
Sik = Cik — oS ik; (5) 
where 
O ff «3 &k, 
bin = “i 
[se gm &. 
The stress-strain relations of an isotropic, incompressible elastic material can then be 
written in the form 
Sik = 2Geix, (6) 


where G denotes the modulus of rigidity. 
In view of (5), the equilibrium condition (2) yields 


Sik,k +O 40% = Sik,zk $5 = O. (7) 


But, according to (6) and (1), 
Sik,k = 2Geix,n = G(Ui,nk + Unix) = GUi,ne, (8) 


since for an incompressible material u,.,.=0 and, consequently, ux,i=Uk,x4i =0. Com- 


paring (7) and (8), we find 
‘~~? = Gui, kk (9) 


Hence 
04 = — Guinn = 9, (10) 
on account of the incompressibility of the material. 
According to (5), 
ik, = Siku tb O,u0ik = Siku, 
on account of (10). Making use of (6), (1) and (9), we transform this in the following 


manner: 
Cin = 2G, = Guin + Ux) = — 20, 


Thus 
nx, + 20, = 0. (11) 
In the case of an incompressible elastic body in equilibrium the boundary condi- 
tions may be given in the form of three functions f;(x) which define the components 
of the forces (per unit area) applied to the surface of the body. The forces f; must, 
of course, be in equilibrium, i.e., the surface integral of f;(x) must vanish for 7=1, 2, 3. 
If , denotes the unit vector directed along the exterior normal of the surface of the 
body, the stress components at the surface must then satisfy the conditions 


CiKnk, = ti (12) 
The values of the surface stresses, in conjunction with Eqs. (2) and (11), define the 
stress distribution in the body and, consequently, also the strain distribution and, to 
within a rigid body displacement of the entire body, the displacement components. 
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On the other hand, the displacement components may be given on the surface 
of the body. These given surface displacements must, of course, be compatible with 
the assumed incompressibility of the material, i.e., the surface integral of the normal 
displacement component uj; must vanish. Elimination of o from (9) furnishes 
Ui,kkt— U1,eei =O, Or, after a change of subscripts, 


Miku — Uri = 0. (12a) 


Eqs. (12) in conjunction with the condition of incompressibility, u;,,=0, and the given 
surface values determine the displacement components. 

2. Stress-strain relations of visco-elastic materials. Equations similar to (11) and 
(12a) may be derived for visco-elastic materials characterized by linear relations be- 
tween the components of stress, strain and their derivatives with respect to time. 

In the case of an incompressible material of the type considered by Voigt? we have 
the stress-strain relations 

Sik = 2Geix, + 2péix, (13) 
where yp is the coefficient of viscosity. 

In the case of an incompressible material of the Maxwell type, we have 

é Sin + (14) 
Ci = ~— Sih a a 
2G 2Gr 
where dots denote differentiation with respect to time, and 7 is the relaxation time.* 

Generalizing, we may consider incompressible materials characterized by stress- 

strain relations of the form 


@* "hadi 0” ite 
(. mee Binh nt a es as) Sik = (0, + 6..3—— F *+* + be (15) 
at™ arm at” ot" 








where @m—1, * * * » @o, Bn, Bn—1, - * * » bo are constants characteristic of the material. 

For such materials two types of boundary value problems may be considered. In 
the first case the surface forces f;(x, ¢) are given as functions of the position x and the 
time ¢; for £=0 these surface forces and their m—1 first derivatives are supposed to 
vanish as well as all stress components and their derivatives up to the order m—1. 
Moreover, at any given time the forces f; must be in equilibrium. If, for #20, the 
forces are analytic functions of time, this implies that the surface integral of any 
derivative 0°f;/dt? must vanish for, say, t=0. The first boundary value problem calls 
for the determination of the stress distribution o(x, ¢) fulfilling these boundary con- 
ditions and initial conditions. 

In the second case the surface displacements u;(x, #) are given as functions of the 
position x and the time ¢; for £=0 these surface displacements and their m—1 first 
derivatives are supposed to vanish, as well as the displacements in the interior of the 


2 W. Voigt, Abh. Gottingen Ges. Wiss. 36 (1899), 47 pp. 

? R. Simha has recently used the stress-strain relations which are obtained from Eqs. (14) by sub- 
stituting the stress tensor ox for its deviatoric part s;x [J. Appl. Phys. 13, 201 (1942) ]. Such stress-strain 
relations imply that, at constant strain, the stress decays exponentially with a relaxation time r which is 
independent of the geometrical nature of the stress. This treatment ignores the fact that viscous flow, 
which is the cause of relaxation, is a response to shearing stresses only. In an incompressible material a 
uniform hydrostatic pressure does not produce viscous flow, and, hence, does not tend to relax. Contrary 
to Simha’s stress-strain relations, our Eqs. (14) reflect this behavior, 
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body and their derivatives up to the order n— 1. Moreover, on account of the assumed 
incompressibility of the material, the surface integral of the normal displacement 
component ui; must vanish for any time. If, for 20, the displacements are analytic 
functions of time, this means that the surface integral of all expressions of the form 
(0°u;/0t?)n; must vanish for, say, t=0. The second boundary value problem calls for 
the determination of the displacements w;(x, ¢) in the interior of the body, fulfilling 
these boundary conditions and initial conditions. 
Let us rewrite the stress-strain relation (15) in the form 


Ps ix = 20¢:x, (16) 


where P and @ denote the linear differential operators 














o™ ont 0 
P= —-+ Gm-1 Pies > ey — + ee, 
ot™ oi™—1 ot 
(16a) 
o” gn-l Ps] 
Q=), + Dn-1 +--+ +b, —+ be 
ot” on Ot 


Starting from the stress-strain relations (16) and repeating the various steps which 

led to the Eqs. (11) and (12a), we obtain 
P(oix,u T 20 ik) = 0 (17) 

and 
O(uin0 — Uxin) = 0 (18) 
as the equations governing the solution of the first and second boundary value prob- 
lem, respectively. 

For example, consider the first boundary value problem for an incompressible ma- 
terial of the Voigt type. Comparing Eqs. (13), (16) and (16a), we see that for this 


material 
re) 


P =i, QO = u—+G. 
Ot 

Eq. (17) consequently takes the same form as for an incompressible elastic material 
(see Eq. (11)). This means that, im the case of the first boundary value problem, the 
stress distribution in an incompressible material of the Voigt type is identical with that 
in an incompressible elastic material under the same instantaneous surface forces. This 
stress distribution does not depend on the past stressing history, although, of course, 
the displacements do. 

This result is readily extended to the case of an incompressible visco-elastic ma- 
terial characterized by a stress-strain relation (16). Consider, for instance, the first 
boundary value problem for a given set of surface forces f;(x, 4) which, in addition to 
fulfilling the conditions stipulated above, are supposed to be analytic functions of 
time for t20. If (x, #) denotes the static‘ stress distribution in an incompressible 
elastic body of the same shape which is subjected to the surface forces f;(x, t), the 
required stress distribution in the visco-elastic body is given by 


oix(x, t) = oix(x, 2). 


4 The term “static” is used here to indicate that, though the stresses ¢;, depend on tas do the forces f;, 
no inertia effects should be taken into account in computing these stresses. In fact, as far as this elastic 
body is concerned, ¢ plays the role of a parameter which need by no means be identified with the time. 
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Indeed, by definition, the stresses ¢ satisfy the conditions (2), (11) and (12) for any 
value of ¢. Since, like the surface forces, these stresses are analytic functions of time, 
this means that they also satisfy the condition (17). The result formulated above for 
the first boundary value problem of an incompressible material of the Voigt type applies, 
therefore, to any visco-elastic material characterized by stress-strain relations of the form 
(16). 

A similar result is obtained in the case of the second boundary value problem for 
an incompressible visco-elastic material obeying stress-strain relations of the form 
(16), if the prescribed surface discplacements u;,(x, ¢) fulfill the conditions formulated 
above and, in addition, are analytic functions of time for +20. The displacements 
u(x, t) then equal the static displacements #;(x, t) of an incompressible elastic body 
of the same shape, subjected to the given surface displacements u;(x, ¢). 

3. Determination of the displacements in the first boundary value problem of 
visco-elasticity. Let us first consider the particularly simple case, where the given 
surface forces can be factored into the form: 


fi = filx)g(d. (19) 


According to what has been said above, the stress distribution which these surface 
forces produce in the visco-elastic body has then the form 


oix(x, t) = Gix(x)g(d), (20) 


where 6;;(x) denotes the stresses which the surface forces f;(x) produce in an incom- 
pressible elastic body of the same shape. Introducing the stresses (20) into the stress- 
strain relation (16), we see that the strains in the visco-elastic body can be written 
in the form 
ex.(x, t) = éx(x)h(d), (21) 
where h(t) satisfies the differential equation 


Oh = Pz, (22) 


while # and its derivatives up to the order »—1 vanish for f=0. As regards the quan- 
tities é;,(x), they are related to the stresses ¢i,(x) by 
(23) 


Sik == Ex, 


where §, denotes the deviatoric part of the stress tensor ¢;,. In other terms, the quan- 
tities €, are the strains in an incompressible elastic body of the same shape and of 
unit modulus of rigidity, which is subjected to the surface forces f;(x). We shall call 
these strains the equivalent elastic strains. In order to obtain the function h(t), all we 
have to do is to consider the response of the visco-elastic material under consideration 
to a simple shearing stress s varying according to s =2g(t). The shearing strain pro- 
duced by this stress equals h(t). The strains produced in the visco-elastic body by the 
surface forces f;(x) g(t) are then obtained by multiplying the equivalent elastic strains 
by the response function h(t). 

Since the differential equations for stresses and strains are linear, solutions of this 
type may be superimposed on each other. Let us, now, assume that our result holds 
good even if, contrary to the assumption made above, the surface forces are not ana- 
lytical functions of time for =0. In particular consider the case when f;=fi(x)g(é, 4), 
where g(£, ¢) is Heaviside’s unit step function defined by 
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0, if ¢< §, 
= &. 


(i, ) = . itt 


Let h(£, #) denote the response of the visco-elastic material under consideration to a 
simple shearing stress s = 2g(é, ¢). Since the surface forces f;(x, t) can be represented 


in the form 


fds,0 f fils, Delt, db, (24) 
0 


the following formal integral representation of the strains produced by these surface 
forces in the visco-elastic body suggests itself: 


wie f guna, &h(E, Odk, (25) 
0 


where é(x, £)dé are the equivalent elastic strains corresponding to the surface forces 
fi(x, £)dé. It can be shown that (25) indeed furnishes the strains of the visco-elastic 
body whenever the surface forces can be represented in the form (24). Moreover, to 
within a rigid body displacement the displacements of the visco-elastic body are 
given by 


u(x,t) = f a(x, E)h(E, tdé, (26) 
0 


where 4;(x, &)d& are equivalent elastic displacements produced in an incompressible 
elastic body of the same shape and of unit modulus of rigidity, by the surface forces 
Fix, £)dé. 

Let us consider the following example: A thin cantilever beam of length Z and 
cross sectional moment of inertia J is clamped rigidly at the end x =0. The beam con- 
sists of an incompressible visco-elastic material of the Voigt type (stress-strain rela- 
tions (13)), and is subjected to the transverse load 


x 
idermalem Se 
2, =C¢ ; 


per unit of length, c being a constant. At first sight, it may seem that the problem 
of determining the bending moments and transverse displacements of the beam is 
outside the scope of our theory, since at the clamped end we have prescribed de- 
formations rather than prescribed forces. However, the system being statically de- 
terminate, the transverse reaction and the bending moment at the clamped end are 
completely determined by the given loads. Consequently, the problem may be con- 
sidered as a first boundary value problem, if we make the usual assumption that 
the distribution of stresses over the end section is irrelevant as long as it leads to the 
resultant and the resultant moment required by the equilibrium of the beam. The 
displacements of an incompressible elastic cantilever beam of unit modulus of rigidity, 


loaded by f(x) =c(1—x/L), are 


cx? 
a(x) = ——— (10L* — 10L?x + 5Lx? — x‘), 
3607L 
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where account has already been taken of the fact that the Young’s modulus equals 3G 
for an incompressible elastic material. 

Now, in accordance with (13), the response h(#) of Voigt’s material to a simple 
shearing stress s = 2? is found from 


22% = 2Gh+2uh, h(0) = 0. 


One obtains 
Mm 


ed a Gt) 
—~Gtlu 
h(t) : a t+ 2 = [1—e |. 


i 


The deflection of the visco-elastic beam is, therefore, given by 


3 2 

u(x, 1) = ———— [10L* — 101% + SL2* — 2°]. E 3 —§+ 8 h< ena). 

360GIL G G? 
The statically determinate bending moments are completely determined by the given 
loads. 

4. Determination of the stresses in the second boundary value problem of visco- 
elasticity. A similar procedure leads to the determination of the stresses in the second 
boundary value problem of visco-elasticity. Consider first the case when the given 
surface displacements can be factored into the form u;=u;(x)g(t), and denote by 
Gix(x) the equivalent elastic stresses, i.e., the static stresses set up in an incompressible 
elastic body by the surface displacements u;(x). Furthermore, determine the response 
function h(t), i.e., half the shearing stress produced in the visco-elastic material under 
consideration by a simple shearing strain g(t). The required stress distribution in the 
visco-elastic body is then given by oxx(x, t) =o;(x)h(t). 

In the general case, the stresses in the second boundary value problem may be 
represented in the form 


on(x, t) = f Gix(x, £) h(E, Adé, (27) 


where &i(x, £)dé are the equivalent elastic stresses corresponding to the surface dis- 
placements w;(x, £)d&, and 2h(§, t) is the response of the visco-elastic material to a 
simple shearing strain 


0, if t<é, 
1, if #26 


g(é, t) = \ 


5. Summary. The solution of the first and second boundary value problems of 
visco-elasticity is reduced to the solution of equivalent boundary value problems of 
elasticity, and the determination of the response of the visco-elastic material under 
consideration to a simple shearing stress or a simple shearing strain. It remains to be 
seen in how far the technique developed here can be applied to the solution of the 
third (mixed) boundary value problem where the surface forces are prescribed on 
part of the surface of the body, and the surface displacement on the rest of this 
surface. 
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THE INTRINSIC THEORY OF THIN SHELLS AND PLATES 
PART III.—APPLICATION TO THIN SHELLS* 


BY 
s 
WEI-ZANG CHIEN 
Department of Applied Mathematics, University of Toronto 


10. Definitions and method of approximation. The method of approximation used 
below is essentially the same as in the case of thin plate theory. We define € to be 
the average reduced thickness of a shell. (We may recall that the reduced thickness 
of a shell is the ratio of its thickness to a selected lateral dimension of its middle sur- 
face). Then for a thin shell, € is a small quantity. This definition of a thin shell is in 
agreement with that of a thin plate given in Part II. 

A thin shell is said to have finite curvature when the smallest radius of curvature 
of its middle surface and the selected lateral dimension are of the same order of mag- 
nitude. Furthermore, a thin shell is said to have small curvature of order b when the 
ratio of the selected lateral dimension to the smallest radius of curvature of its middle 
surface is of the same order of magnitude as e’, where 21. Thus a thin plate may be 
regarded as a thin shell of small curvature of order ~. 

We consider a family of «! shells of the same material with diminishing reduced 
thickness, each in a state of stress under (i) external forces applied at the edge, 
(ii) surface forces and (iii) uniform body forces. We assign to each shell a value of a 
parameter e (0<¢€<.e,) denoting the average reduced thickness, so that the thickness is 


2h = 2ch(x!', x”). (10.1) 
The quantity € is supposed to be small, but the basic idea of the method is that we 
seek solutions valid for all € in the range 0<e<«. In this theory, € is the only small 
quantity. All quantities occurring (except Poisson’s ratio o) are functions of €. No 
quantity is small unless it tends to zero with e. 


For the greatest generality suppose all quantities to be power series in e. Thus, 
supposing the middle surface itself to depend on e, we have 
io) i) 
aap = 2. A (8) ape’, bas = =. Bs) asé’, (10. 2a, b) 
s=b 


&=U 
where 3 is either zero or a positive integer. a;s)ag and bys)as are functions of x*, inde- 
pendent of e. For )=0, we are dealing with thin shells of finite curvature, while for 
b=1 we are dealing with thin shells of small curvature of order b. 
. : rt wm mn 4 ; vere ea 
Furthermore, we shall represent Q', P‘, Xjq, T°, T*°, L*°, pas, das by power series 


as in Part II; 
= > Oe, Q7 = > Oz", (10. 3a) 


* Received June 12, 1943. Parts I and II of this paper appeared in this Quarterly, 1, 297-327 (1944), 
and 2, 43-59 (1944). 
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po = 3 Phe, Pa = Pees, (10. 3b) 
&=no s=k 
x, = £ Xoo, X% = © Xone (10. 3c) 
po = T° Ther, ie = ive, To = Te, (10.40, b, ©) 
s=l s=u s=l 
Pas = > Ps) ap€’; Gas = = F (s) as€". (10. 5a, b) 


Here k, ko, ”, "0, j, jo, t, u, l, p are integers greater than zero, and q is zero or a posi- 
tive integer. The case g =0 corresponds to problems of finite deflection. The quantities 


Q°.), 0%, Pt ete. are functions of x*, independent of e. 





p-values 
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Fic. 4. Classification of problems of thin shells with finite curvature (b =0). 
p=order of extension of middle surface. 
q=order of change of curvature of middle surface. 
b=order of initial curvature of middle surface. 


Then the problems of thin shells can be classified by assigning integral values to 
b, g and b. With 9, q, b given, the values of ko, k, mo, ”, jo, 7 in (10.3) are fixed by the 
condition that X?;,yo, XGyo Pes» Pe» Qf)» Qf% should contribute to the principal 
parts of (6.34), (6.35), without dominating these equations to the exclusion of Pag 
and qas. The values of t, u, 1 of 78, L*8, T° are immediately fixed through the ex- 
pressions (6.29), (6.30), (6.31). With p, g, b, k, ko, j, jo, m, mo fixed, the equations of 
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equilibrium and compatibility in the first approximation are immediately obtained 
by substituting (10.1)—(10.5) into (6.34), (6.35), (6.43), (6.44), and picking out the 
principal terms in ¢ from the resulting equations. This gives us six differential equa- 
tions in six unknowns P;»)as and Qig)ag. For the various combinations of values of 
p, g, b, the forms of these differential equations fall into several types. The classifica- 
tion of these types will be given below. 

11. Classification of all thin shell problems. The classification of the problems of 
thin shells with finite curvature (b =0). The following is a complete classification of the 
problems of thin shells with finite curvature (b=0) based upon assigned values of p, q. 
The classification is shown graphically in Fig. 4. 

It is found that the (p, g)-points in the diagram (¢20, p21) are broken up into 
eight groups by the division lines AB, OC and the p-axis. For g=0, the principal part 
of (6.34) or (6.35) takes three different forms depending on the position of the point 
on the p-axis relative to the point A, while the principal parts of (6.43) and (6.44) 
are the same for all values of p. For g21, the principal part of (6.34) or (6.35) takes 
three different forms depending on the position of the (p, g)-point relative to the line 
AB, and that of (6.43) or (6.44) takes three different forms depending on the position 
of the (p, g)-point relative to the line OC; each of these forms is different from that 
for g=0. It follows that the (p, g)-points are divided into eight groups and so the 
complete classification of all problems of thin shells of finite curvature involves con- 
sideration of eight types (Types SF1—SF8). (The letter ‘S’ denotes shell, while ‘F’ 
denotes finite curvature.) 

In order to save space, we shall not discuss these types in detail. The results for 
these types are summarized together with those for thin shells with small curvature 
in the tables in the Appendices. The principal parts of the equations of equilibrium 
and compatibility are shown in Table III, and orders of magnitude of the external 
forces and the principal parts of the macroscopic tensors in Table IV. 

The classification of the problems of thin shells with small curvature (b2=1). The 
following is a complete classification of the problems of thin shells with small curva- 
ture based upon the assigned values of b, p, g. The classification is shown graphically 
in Fig. 5 (for b=4), Fig. 6 (for b=2), Fig. 7 (for b=1). The case b=4 is typical of the 
cases 3S50< ~. 

We shall now explain Fig. 5. We see that the (, g)-points are broken up into 
27 groups by the division lines and the p-axis. Of these division lines, the line B’BB’’ 
(i.e.,g=b=4) is the most important. It divides the (p, g)-plane into three main re- 
gions. For any point on B’BB”’’, the curvature in the unstrained state and the change 
of curvature during the strain are of the same order of magnitude (¢=) =4). For any 
point on the left of B’BB’’, the magnitude of the curvature in the unstrained state 
is smaller than the magnitude of the change of curvature (¢<b=4), while for any 
point on the right of B’B’’, the magnitude of the curvature in the unstrained state 
is greater than the magnitude of change of curvature (¢>b=4). 

For g=0 (i.e., on the p-axis) in Fig. 5, the principal parts of (6.34), (6.35) take three 
different forms depending on the position of the points on the p-axis relative to the 
point A, while the principal parts of (6.44), (6.43) are the same for all points on the 
p-axis. For 1q<b=4 (i.e., in the region between the p-axis and B’BB’’), the prin- 
cipal parts of (6.34) or (6.35) or (6.44) take three different forms depending on the 
position of the (, g)-point relative to the division line AC or AB or OD respectively, 


1944] INTRINSIC THEORY OF SHELLS AND PLATES 123 
while the principal part of (6.43) is the same for all the (, g)-points in this region. 


It follows that the (p, g)-points in the region on the left-hand side of B’B’’ are divided 
into 11 groups (Types SS1—SS11). (The letters ‘SS’ denote the shell with small curva- 


ture.) 
p-values Bt 
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Fic. 5. Classification of problems of thin shells with small curvature (b =4). 
p=order of extension of middle surface. 
q=order of change of curvature of middle surface. 
b=order of initial curvature of middle surface. 


For g=b=4 (i.e., on B’B"’), the principal parts of (6.34) or (6.35) or (6.44) take 
three different forms depending on the position of the (, g)-point relative to C or 
B or D respectively, while the principal part of (6.43) is the same for all points on this 
line. Furthermore, for g>b=4 (i.e., the region to the right of B’B’’), the principal 
parts of (6.34) or (6.35) or (6.43) or (6.44) take three different forms depending on 
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the position of the (/, g)-point relative to the division line CG or BE or B’H or DF 
respectively. It follows that the (p, g)-points on the right-hand side’ of B’B’’ are di- 
vided into 9 groups (Types SS19—SS26, SS10). It should be noted that, as far as the 
principal parts of (6.34), (6.35), (6.43), (6.44) are concerned, the (p, g)-points lying 
between the lines DF and JCG are regarded as one group (Type SS10). Therefore, 


p-values Bt E F 
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q-values 


Fic, 6. Classification of problems of thin shells with small curvature (b =2). 
p=order of extension of middle surface. 
qg=order of change of curvature of middle surface. 
hb =order of initial curvature of middle surface. 


together with the groups on the left-hand side of B’B’’, we have in all 25 groups of 
(p, g)-points in Fig. 5. And consequently the complete classification of all problems 
of thin shells with small curvature of order }=4 involves consideration of 25 types 
(Types SS1-SS11, SS13-—SS26). 

The general appearance of the classification diagrams for any b satisfying 3b < « 
is the same as for b=4. An increase of b makes the line B’B”’ shift to the right, while 
a decrease of b makes it shift to the left. On examining the various groups of (p, q)- 
points in these diagrams (for any integral value of 6 in the range of 3Sb< ~), it is 
found that the corresponding groups occupying the same relative positions with re- 
spect to the division lines possess the same set of equations of equilibrium and com- 
patibility in the first approximation, and so belong to the same type of problem. 
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Therefore the complete classification of all problems of thin shells with small curva- 
ture of order 3b< & involves consideration of 25 types only. 

For b=2 (Fig. 6), the situation is almost the same as in Fig. 5, but with the groups 
SS9, SS11 missing. The other groups are the same as those shown in Fig. 5 for )=4, 
and so no extra types arise. 

For b =1 (Fig. 7), the situation is only slightly different from those in Figs. 5 and 6. 
Instead of the two separate division lines JDF and ICG for Eqs. (6.34) and (6.43) in 
Figs. 5 and 6, we have one common division line D’F’ for both equations. Further- 
more, the triangle formed by the division lines JD, DC, IC in Figs. 5, 6 collapses into 


p-values Bre E Ff 
9 

















0 1 2 3 4 5 6 7 8 
ee 


q-values 


Fic. 7. Classification of problems of thin shells with small curvature (b=1). 
p=order of extension of middle surface. 
q=order of change of curvature of middle surface. 
b=order of initial curvature of middle surface. 


an isolated point D’ in Fig. 7. Thus instead of 25 groups in Fig. 5, or 23 groups in 
Fig. 6, we have only 15 different groups. Among these groups, 13 belong to the types 
already mentioned in the case 3Sb< © (Types SS1-SS3, SS13, SS16-SS21, SS24- 
SS26); the other two are Types SS12, SS27. 

On comparing the classification of (p, g)-points on the left-hand side of B’B’’ in 
Figs. 5, 6, 7 with that in the corresponding region of Fig. 3, it is found that they 
are identical with each other. In fact, for these types, the equations of equilibrium 
and compatibility in the first approximation are identical with those stated in Table I 
(Part II) for the corresponding types of thin plate problems. Therefore, we have the 
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following important conclusion: A problem of a thin shell with small curvature of order b 
is effectively equivalent to a problem of a thin plate in the first approximation, if q<b, 
i.e., if the change of curvature is greater than the curvature of the shell in the un- 
strained state. 

It should be noted that for b= ©, Fig. 5 becomes exactly Fig. 3 for the thin plate 
problem. 

The results are summed up as follows: 

(i) The complete classification of the problems of thin shells with small curvature 
of order b21 involves the consideration of 27 types (Types SS1—SS27). 

(ii) Among these 27 types, 11 are equivalent to problems of thin plates; the char- 
acteristic of these types is g <b. 

(iii) When 6=1, these are two types (Types S$S12, SS27) of particular interest. 

We shall not discuss all these types in detail. The discussion of Type SS12 will 
serve as an example. The results for all types are summarized in tables in the Ap- 
pendices. The principal parts of the equations of equilibrium and compatibility are 
shown in Table III, and the orders of magnitude of the external forces and the prin- 
cipal parts of the macroscopic tensors in Table IV. 

Before entering on the detailed discussion of Type SS12, a useful result for small 
curvature (b21) will be mentioned. On substituting aas, bag from (10.2a, b) into 
(6.39b), it is found that the lowest power of ¢€ in the resulting expression is e®. The 
corresponding coefficient gives rise to the equation 


Ro) paby = 3(8(0)py,a8 + 8(0)a8,9y — 2(0)p6,ay — 8(0)a7,8p) 
+ af) { [ov, tJeolaB, 5Joo — [p8, tlaglay, 5Jo} = 0, (11.1) 


where the Christoffel symbols are calculated for ayo)as. Eq. (11.1) expresses the fact 
that in the case of small curvature, the curvature tensor vanishes in the first approxi- 
mation. Hence the order of the operations of covariant differentiation with respect 
tO @o)ag is immaterial; this result will be found very useful later. 

12. Detailed discussion of type SS12 (b,=q=1, p=2) and its applications. General 
equations. By the condition that in the first approximation, (6.34), (6.35) receive sig- 
nificant contributions from P%.), P&, Xgooy XGyoy Ole)» Qa», We must have 


no = 4, jo = 3, ko = 2, 
n= 3, j = 2, hoo 3. 
By substituting the € series into (6.34), (6.35), (6.43), (6.44), it is found that the 


lowest powers in € occurring in the resulting equations are respectively e*, e*, €', €’. 
The corresponding coefficients give rise to the following equations: 


(12.1) 


pyr i. yar 7. 2 yr 7. 
— AGI B(1)9r1P 2) ak — 2A) G1)9rxP 2) ak + FAG (ay ea’) tor 
ao 


i —— : 2(1 — 2c) — 
+ Pu) + 2X (3)[0] h+ (Qsyh) at “1 = HQah 
ao 





a 

1—2¢ ‘aia 
+ 1 ew 21) rA(Q)eayh = 0, (12. 2a) 
— — 


- » og o 
2Afsi) (Peyerh) 1p + P& + 2X Gyoih + — al (Qh) in = 0, (12. 2b) 
ao all ao 
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nasi = 0, (12. 2c) 
a0 
nisin P arerias + mneyaarer@aras + (bh — 4H wad)aaras = 0, (12.2d) 


ali) 


where ap under stroke indicates covariant differentiation with respect to the tensor 
A(o)ag and x*. The other symbols represent 


1 

Ais = ——— (oatia + (1 — o)aiialo), (1h. 

= a? 
nip) = €78( ao) 1/2 an = det. (Aco) sa); ell = ¢72 = 0, 62 = — ¢21 = 1, (12. 3b) 
bi) = avaobana, Ha = taba.’ (12. 3c) 

The macroscopic tensors (6.29)—(6.31) can be written as 
Te? = { 24% Pex wah Mii Fe = a Qh ite + O(¢), (12. 4a) 
L? = sic las ale + O(€), (12. 4b) 
To? = {3A (Qansh*) ix + Oyh}et + O(E). (12. 4c) 
bal’) 


Equations (12.2a, b, c, d) are six equations for the six unknowns py), and q,1 sa. 

Since by (11.1) the order of the operations of covariant differentiation with re- 
spect to @c@) a is immaterial, (12.2c) implies the existence of wa) such that 

Q(1)ap = Wc1)jas- (12.5) 
40 

Thus the determination of qa is reduced to the determination of the single func- 
tion Wa). Furthermore, instead of using p,2)ag as the rest of the unknowns, we may 
use 7#. By definition, Tée* is the principal part of the macroscopic tensor T%, 
namely, by (12.4a), 


TH = 2ACi. 'P (2) ah a al 2h. (12.6) 
This is a symmetrical tensor; so it has only three independent components. Substitut- 
ing (12.5), (12.6) into (12.2a, b, d), we have 
— TEwayia — 4bayaT® + FAC (Wain) py + Plo 
ao ao 


+ 2X ayo) + (Qik) ix + 2H Qh + wan nao = 0, -  (12.7a) 
Tie + Pi + 2X tot = 0, (12. 7b) 

nfo { (1 + 0) aco) x80) 78 — TA (0)p78(0) x8} is TE i + 240) Qta ian 
+ moi «1 189 Lap + (bf — 4H aio) co Job = 0, (12. 7c) 


Equations (12.7a, b, c) form a set of four equations for the four unknowns wa) 
and 7%. 
Special case. The following special case is interesting. If 
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Pi) = Xe = 9, (12.8) 
then by (12.7b) there exists a stress function x,3) such that 
an. Br 
TY) = No N{O)X (3) | a (12.9) 
ao 


Here x,3) is a function of x*, having properties similar to those of the Airy function 
in the thin plate theory. Substituting (12.8), (12.9) into (12.7a, c), we have 


1A. dp / 2 APryTr/, i 
= 31/9) 1/0) ( 2 (1) | rs + BD, 78) X02) hp + 3A (W (1) xrh*) |p 


ao ao ao 


+ Py + 2X Gok + 2H Qh + alwayinQah = 0, (12. 10a) 
40 


B » Ji 
{ catpalo) — (1+ oa (0) 2( * VG X (3) ati | xp + wate) Oe) | xr 
U ao 


a TH 
” +r X 
+ nfo nfo (1) loa i (biG) ~~ 4H (1)a(o))W 1) | wa a= (). (72. 10b) 
ao a ao 


Equations (13.10a, b) are two equations for the two unknowns x) and wa). These 
equations are valid in general for a shell of non-uniform thickness. For the case of 
uniform thickness, (12.10a, b) are immediately simplified to the forms 


™r _ dp | wy rd 
— 3n10 Ni (2W (1)\ x8 ao Bi) x8) X (3) [dp + Dajiy 4(o)W (1) | rydd 


A FY ay 
0 70 i, 07 ™r Je 4 
a Pi) oo 2X (3)[0/4 a 2H, Vo h a aw W (1) | rxQiayh = 0, (iz. 11a) 
ao 
y rd yom 7 dy 
€(0)8(0)X 3)|ryvA8 cha QO 2) |r + hnipno Wi pyW (1) | ra 
ao ao ag ay 
E » 

a h(4H yay — bé))Wain = 0, (12.11b) 


40 
where D is the reduced flexural rigidity, as given in (9.14). Applications of these two 
equations will be discussed below. 
A circular cylindrical thin shell with small curvature and uniform thickness under 
end thrust and normal pressure. We shall assume that the external forces and the edge 
loading are such that the problem is of Type S512. Furthermore let us assume that 


rfayo) = Qo) = 0. (12.12) 


We have in mind the case where body force is negligible and where the shell is loaded 
normally on one side only. A number of terms disappear from the equations of equi- 
librium and compatibility (12.11a, b) for Type SS12. Thus if we write these equations 
in terms of the small principal parts instead of in terms of the finite coefficients of the 
lowest power in €, we have 


ng Anji; (2073 + bss)x - + Da*™aw ryrd + a _ 0, (12. 13a) 
ata y jeans + Anginigw)..wixs + A(4Ha” — b™)win = 0. = (12. 13b) 
a a a a 


Here a under a stroke indicates covariant differentiation with respect to the tensor aa 


and x; also 

2h? 
ae 4H = a.sb”. (12.14) 
~ 34 — o’) 
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Let us choose the set of intrinsic rectangular Cartesian coordinates on the middle 
surface so that x =x! is the distance measured along the generators of the cylinder and 
y =x’ is the distance measured perpendicular to the generators. Then we have 


a), = a" = an» = a” = I, ay: = a’? = 0, 


(12.15) 
by = b'! = bis = b'? = 0, boo = b” = 2/R, 


where R is the radius of curvature of the cylindrical middle surface. In these coordi- 


nates, Eqs. (12.13a, b) become 


1 
DAAw + (2W, 2X, 2y — Wi22X.vv — WiyyX,22) — phat ee, (12. 16a) 


) 1 
AAx + 2h(W,22Woyy — W.ryW, ry) + 2h Pie = (0, (12. 16b) 


where subscripts preceded by a comma denote partial differentiation. If we let R 
tend to infinity, we get the von Karman equations for a flat plate. The equation 
(12.16b) was recently obtained by von K4rman and Tsien [1] in their treatment of 
buckling of a thin-walled circular cylindrical shell under compression on the two ends. 
If we apply the operators AA to (12.16a) and (1/R)0?/0x* to (12.16b) and add the re- 
sulting equations, we obtain 

2h 2h 
DAAAAWw + — W.erzz + — (W.22Wiyy — WzyW zy) ,22 


= AA(P® + 2w, syX.2y — WoezX.vy — WiyyX,22)» (12.17) 


This is the equation of equilibrium used by von Karman and Tsien, except that they 
omit the term 

2h 
bie (W, 22W yy — all al.com (12.18) 


This term is important when the deflection is 
comparable with thickness. However, it seems 
simpler to treat the problem directly by means 
of (12.16a, b) instead of using the higher-order 
equation (12.17). Equation (12.16a) appears 
to be new. 

A small segment of a thin spherical shell un- 
der external pressure. We shall assume that the 
solid angle of the segment is small, so that the 
curvature is small; we shall assume it to be of 
the same order as the thickness, so that }=1 
(cf. section 10). We shall use spherical polar 
coordinates as in Fig. 8, so that on the middle 
surface in the unstrained state we have 


ds* = R°dé? + R? sin? @dg*. (12.19) 








Since @ is small, we write Fic. 8. 
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ds? = R*d0? + R*6*de?. (12.20) 
If we put 


xi=§, x=, (12.21) 


the components of the first and second fundamental tensors are given by 


ai = R*, as = R6*, aw=0, at =1/R’*, a” = 1/R?, a? = 0, (12.22) 
bi, = 2R, be = 2R6?, b= 2/R*, b?? = 2/RO?, bi = Bb? = 0. (12.23) 
Futhermore, we have from (12.22), (12.23) 

H=1/R, a= R*. (12.24) 


All the Christoffel symbols are equal to zero, except 


4 2 ; 
: } wa’, {, } = 1/0. (12.25) 


We shall assume that the problem is of Type SS12. Substituting (12.21)—(12.25) into 


(12.13a, b), we have 


1 ( ‘a 1 1 

on ed Ws + Ge + ted + 2( wo —-— we) (x1 i 7**) 
R402 | p 6 

. ? 


+ (Wwe + RO? + ow s)x.at + DAAw + P®°=0, (12.26a) 


2h 1 . 
AAx + aes | Wm + bw) — (0 mee we) \ 
R‘9? i] 


1 1 \ 
2 a ee ) = 2 
+ 2h ‘ms Wo + — (wee + walt 0. (12. 26b) 
Here A is the Laplace operator 
1 0? 1 0? ; @ : 2 0 1 0? 
= cele (12.27) 





i a a ee 
R? 06? R*0? dg? R*0 00 =§=R*6 00 «00 R*0? dg? 


Equations (12.26a, b) are two nonlinear partial differential equations for two un- 


knowns x, Ww. 
We suppose that the problem has rotational symmetry. Then w, x are independent 


of gy, and (12.26a, b) reduce to the form 


d did _ <dw 1 d/dw dx Rd dx P°@R* 
j -— )-. (« )+ 0, (12. 28a) 


oe” eee [ads aera Pattee rds ess Bi 
do dé 6 dé dé D dé\d6 dé D dé dé D 


@.a@ 14 .& d { dw\? d dw 
—§— — <0 + no () + mar =) 
dé d@ 6 dé dé d6\, dé dé dé 





II 


0.  (12.28b) 


The equations can be integrated once giving 





1944] INTRINSIC THEORY OF SHELLS AND PLATES 131 


d 1d dw 1 dw dx R_ dx P9?R4 


i cee nc aaa 
dé 6 d@ dé D do dé D dé 2D 
d 1d _ «ax (= 


do 6 do dé d0 





= constant, (12.29a) 


: dw 
) h+ 1a —— = constant. (12. 29b) 


Since dw/dé@ vanishes for 6 =0, the constants are zero. If we introduce the quantities 


dw 1 dx 





1 
me oe eh ak ee 12.30 
‘ R dé R?* dé ’ 
the equations can be further simplified to the form 
d*a 4 da a R? B+ P°R® g2 0 (12. 31a) 
rice auiat denis eae ead sais “= Q, ° a 
de | d0 6 D- 
ap dp 8B 
8 4 —— + h(a? — 6?) =0. (12.31b) 
dé? dé 6 


The quantity a is the slope of the meridian line 
in the strained state (Fig. 9). The significance of 
the quantity @ is that 6/0 is the radial membrane 
stress (tension). Equations (12.31a, b) are the 
fundamental equations for the determination 
of the buckling pressure of a small segment of 
spherical shell. 

If we assume that the first and second terms 
in (12.31b) are negligible in comparison with the 
other terms, then we can solve (12.31b) immedi- 
ately for 8. Substituting the resulting expression 

















for 8 into (12.31a), we have Fic. 9. 
, d?a 4 da a hR? dala? — 62) P92 R8 (12.32) 
— rm — = alas — ‘> — — r » 
de?—s« 6 D 2D 


This is the equation used by von K4rm4n and Tsien [2] in their treatment of buckling 
of spherical shells by external pressure. It should be noted that the neglect of the first 
two terms in (12.31b) is a rough approximation. Actually the first three terms in 
(12.31b) are of the same order of magnitude. 
Furthermore, if we introduce 
r = Ré, (12.33) 


Eqs. (12.28a, b) can be written in the form 


iid d 1 d/f/dwd 1 d d P% 
= << 8,27 (eS) 4 =0, — (12.34a) 
dr dr vr dr dr D dr\dr dr RD dr dr D 


s #€ 3 & @ d fd ae © dw 
anemone n—(=) = <(; —) =o, (12.34b) 
dr dr-vr dr ar dr\ dr R dr dr 
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The quantity ¢ is the radial distance measured along the meridian line from the center 
of the shell. We see that these two equations are the same as the corresponding von 
Karmd4n equations for the circular plate under symmetrical loading [3], with the ex- 
ception of the terms proportional to 1/R; this is evident if we make R infinite in 
(12.34a, b). 

A summary of the whole paper was given at the end of the first section (Part I). 


APPENDICES 


(iii) TABLE III.—Table of the equations of equilibrium and compatibility of thin shell problems. 








(6.34) (6.35) (6.44) (6.43) 
Types b qg p One Ser ee eee samt i cian : (es 
0 0 ) 0 i] 0 0 a a aa a a ) ) 
12 13 18 13 18: 18 WB a8 18 | aS aS as 1S 1S 18 8 8 TS S| Jon Jos 
SS1 } 21 0 1 x . = x : S x x 
SS2 >1 0 2 a oe a ae : 2s 2 x x 
SS3* | 21 0 >2 : es = 2 b 2-S x x 
Soe i 3 1sq<b i x x = x x x x 
SSS {| #2 1 2 zg s x x x x S = x 
SS6* =2 2q+1 | x x x x | x x 
SS7* =2 2g+2 x x zzz 2 x x 
SS8* >2 >2¢+2 | x x £evey x x 
SSO | ES 2 | = =z  - x x x x 
f 

23 <p <2 
SS10¢ 2<p <2q 

>2 2+¢—b<g<q+b x x x x | x x 
SS11 23 2<q<b 2q x x x x te: x 
SS12 1 i 2 g =z zs = = Zi x [= 2 x x 
SS13 >1 b 1 x x x x x|x x i< x 
SS14 2 b 2 |x x cess zis x |x x 
SS15 >2 b 2b x x x x is = x x 
SS16* > 1 b 26+1 x x x x | x x x 
SS17* > 1 b 2b +2 x x : a2.2°8 4 2 x x 
SS18* >1 b >2b+2 x x rea BS x x x 
SS19* >i >b <q-—b x x x x x x 
SS20* + >b q—b x x x |x x x x x 
SS21*| 21 >b qg—b+1 x x x |x = x x 
SS22 =2 >b q—b+2 x x x ei = x x x 
SS23 | 22 >b qg+b x x x x x x x 
SS24* 21 >b qt+b+1 x x x x x x 
SS25* > 1 >b q+b+2 x x x x x =z x x 

| 

SS26*| 21 >b >q+b+2 | x x x . x x| x 
SS27 1 >b q+1 x x x x |x | x x x 
y } 
SF1i | 0 0 1 x eS = a x |x zs 3 x | x x x 
SF2 0 0 2 2s @ 24-&8£&.& 212 2 Ee xO x x x 
SF3* 0 0 >2 | i ee fe | Ss s 2 22 = x x x 
SF4 | 0 =1 q x x siz x x |x x x x x 
SFS | 0 21 q+1 x x x |x x x x x 
SF6 0 21 q+2 3 x x xj|x x z = x 
SF7* 0 >1 >¢+2 } x x x x x =e] x x 
SF8 0 >2 <q } x x ee x ae x x 





In this table, the following notation is used: 
The terms occurring in the first equation of equilibrium (6.34) are 


0 pywwrd 


0 yar 0 2 yar 3 3 
qr; —_— 2Af; dp yP mh, Ts = 3 (i) (qh J lev I — (3) FreQoyQrsh ’ 
a 


1944] INTRINSIC THEORY OF SHELLS AND PLATES 133 








' 1—2¢ 
f= P°+2Xph+ Map, R= ang?Qh, I) = — Ab, pah, 
a - ¢ 
0 oO 0 weeds 0 1 — 20 
7 = — pAb, ybusQreh’, I; = Afé) QrwPpyQrsh’*, = 1 2HQh. 
—o 


The terms occurring in the second and third equations of equilibrium (6.35) are 
If = 2A63(Peh)in TE = Baar AU Qnsh’) ip — ABT (areaesh’) in 
Ig = P* + 2Xph + a Ps a°“(Q°A) ip Tf = (a?qnas + 2a%qr7)Qh, 
Is = AGS (Buadroh!) |p + SAU a" ry (Qnah®) oy Ig = (2Hal + bf)Q%h. 


The terms occurring in the first equation of compatibility (6.44) are 


0 _ By a 8 
Jj — Zn} M0] Po | abs J _ nip) Mio} Dexa» 
a 
a , 
J3 = 2a®pasK, Jt = — (4Ha® — b%*)qas. 


The terms occurring in the second and third equations of compatibility (6.43) are 
Jai = 2n{o}Qas |) J a2 = nij}bs28™(Par|y + Parla — Pay): 
a a a a 


On account of the conditions which hold in the various types of problems, some of 
these terms may be negligible in comparison with others. Table III shows by the 
symbol ‘x’ those terms which are to be retained in the first approximation for the 
various types. (The over-determined problems are denoted by *.) Thus for example, 
for problems of type SS1, we have the following equations of equilibrium and com- 
patibility in the first approximation: 


P+n4+M%=0, w+mt+m=0, JS=0, Ja=0. 


These equations are written in terms of the small principal parts instead of in terms 
of the finite coefficients of the lowest power in e. 

(iv) In Table IV, the following notation is used: 

The terms occurring in the expresion (6.29) for the membrane stress tensor T** 
are denoted by 


| _ 2A" ph, ™? a Se Ag g .qnah', 
o 

ead _ Si ah, TH its Ae™ bysq rh’. 
—_ — 


The terms occurring in the expression (6.30) for the bending moment tensor L*4 
are denoted by 


ag 2.8 anhdd 3 
Ly 31[0}@ rp Aci) geh’, 


8 ‘ at) 
| bs 2nfiar AG ” Dy sPpyh® 


II 


o 4o 
+ ————-a Sa (—*_ xo" — 4Xjq — 203,) = bsort we. 
ti-¢ °° t =~ A 
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TABLE IV.—Table of the external force system and the macroscopic tensors 
for various types of thin shell problems. 











Y ead | aad | T? 0 
Types no n je j ko | k fe i. » a 
2 wt gt et on a owe wees 
S51 2 2 i i 1 1 2 x x 3 :< 2 7 
2 2 i i 1 2 2 x x 3 x | 3 x x 
SS2 3 3 2 2 2 2 3 x x 3 |x ee: x 
SS3 3 3 2 2 2 2 3 x x 3 | x 3 x x 
aT eo eee ee eee 
q@+2 2 q+ qgt+2 2 x x qr: x | @+3 x x 
SS5 4 3 3 2 2 3 3 x x 4 | x | 4 x x 
SS6 | @q+3 g+2|2¢g+1 | 2¢+1 q@+2 | 2¢+2 x x q+3 | x | @+3 s « 
SS7 q+3 2g+2| g+2 | g+2 |2g+2 |2¢+3 | x x x q+3 | x | @+3 |x x 
SSB | a+3 2g+2| g+2 | g+2 |2¢+2 | 2¢4+3 | x x q+3 er lets |x =x 
SS9 g+3 q+2 2 2 q@+2 3 x x qg+3 x q+3 x x 
S$S10 | 9 +3 g+2 p Dd q+2 p+1 x x q+3 x q+3 x x 
SS11 @+3 | 2¢+1 qg+2 2q |! 2¢ | @+2 | 2¢+1 x x q+3 | x | q+3 x 
SS12 4 3 3 2 2 3 3 x x 4 x 4 x x 
ssi3{| +2 2 | b+1) 1 1 | b+1 2 x x b+3 | x | b+2 | x 
‘ 1} 642 2 b+1 i 1 b+2 2 x x b6+3 | x b+3 x x 
SS14 +3 3 b+2 2 2 b+2 3 x x b+3 x b+3 x x 
SS15 | b+3 26+1 b+2 2b 2b | b+2 26+41 x x | 643 x | b+3 x x 
SS16 | b+3 | 2b+2 | b+2)2b+1 [2641 | b+2 [2642 | x x b+3 | x | b+3 | x x 
SS17 | b+3 | 2643 | b+2]2b+2 |2b4+2 | b+2 | 26+3 | x x x | 6+3 | x | b+3 |x x 
$518 | b+3 [2643 | b+2/2b+2 |2b+2 | b+2 | 2b+3 x x x | 64+3 |x b4+3 |x x 
(lb+p+1] pti | b+0| >» p | b+p | pti | x x b+p+3] x |b+p+1| x 
$S19{/b+p+1| pti | b+n| p d |b+p+1) pti | x x b+p+3 x lb+p +2] x 
\/b+p+1 P+1 | b+2 y | #2 lb+p +2 p+1 x x b+p+3 x lo+p +3} x x 
b+p+i] p+1 | b+p) vo» | vb | +p | pti} x x b+p+3| x x [b+p+1) x 
SS20}\b+p+1) p+1 | b+0] » | pb |[b+o+1] pti | x x b+p+3| x x |b+p+2) x 
b+p+i| p+1 b+p p ? lb+-p+2 p+1 x x b+p+3}| x x en . = @ 
ssa Jotett] o+1 | b+] & | & | b+ | ptt | x x b+p+2) x Jb+p+1) x 
C08" "\b+o+1| o+1 | b+0) @ | © |broti] pti | x x b+p+2| x b+p+2]| x x 
SS22 | g4+2 |q—b+3!| ¢+2/¢—b+2\¢—b+2) g+2 |q—b+3| x x q+3 | x | g+3 | x x 
SS23 a+3 |q+b+1| ¢g+2! g+b | a+b q+2 \q+b+1| x x q+3 | x q+3 x x 
SS24 a@t+3 |¢@4+642) ¢+2\¢+b+1/¢+b4+1) g+2 |¢+b4+2!) x x q+3 | x q+3 x x 
SS25 @t+3 j¢4643) ¢4+2\'¢+642'¢4+6b4+2) g+2 |¢+b+4+3) x x x q+3 | x | @+3 x x 
SS26 at3 j|¢@4+b4+3) ¢4+2,¢+64+2'¢+642) ¢@+2 (|¢+b+3 x x q+3 x | g+3 x x 
SS27 | q+3 | g+2 | q+2) g+1 q+1 qt+2 | g+2 | x x qt+3 | x q+3 |x x 
or 2 2 1 a 1 | 2 | x 3 | x | 2 |x 
2 2 1 ; 2 3 z 2 2. y x ae Ss te: 
SF2 3 3 2 a 2 3 . ££ = 3 x s ts << 
SF3 3 3 2 teow Se 3 s+. 3 x . ie 2 
( q+1 q+1 q qg qd qg g+1 x x | @+3 x q+1 x 
SF4 q+1 q+1 q qg q q+1 q+! x x q+3 Ss « q+2 x 
q+1 g+i1 q qg q q+2 q+1 x x q+3 | x x qa+3 x x x 
SFS f q+2 q+2 q+i1| g+1 q+1 q+1 q+2 x x q+3 x | q+2 x 
|} @+2 qgq+2 gq+i} qg+i q+1 q+2 q+2 x x q+3 x q+3 x x 
SF6 q@+3 g+3 qg+2) q+2 g+2 q+2 q+3 x x x qa+3 x q+3 x x 
SF7 q+: q+3 q+2| g+2 qg+2 q+2 q+3 x x q+3 x q+3 x x 
(| p+1 | p+1 | b p p p p+1 | x x p+3 x | pti | x 
SF8 4| p+1 p+1 p b p p+1 p+1 x x p+3 x | p+2 x 
\| o+1 p+1 Dp p i p+2 p+1 x x p+3 x p+3 x x 








The terms occurring in the expression (6.31) for the shearing stress tensor T° are 
denoted by 
a a ward 
Ty = Qch, Tz = $ Ath (Qrsh*) |, 
a 
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T3 = 2A%3"(byapprh®)|2 + F(4HP* + BSP*)h? + $HX ght 


Further mor e, 


ny =order of sum of the normal forces acting on the upper and lower boundary 

surfaces, or order of P®, 

n =order of sum of the tangential forces acting on the upper and lower boundary 

surfaces, or order or P*, 

jo =order of normal component of body force, or order of X{Q, 

j =order of tangential component of body force, or order of X{qj, 

ko =order of difference of normal forces acting on the upper and lower surfaces, 

or order of Q°, 

k =order of difference of tangential components of forces acting on the upper and 

lower boundary surfaces, or order of Q+, 

t =order of membrane stress tensor JT, 

u =order of bending moment tensor L*%, 

1 =order of shearing stress tensor T°. 

This table gives (a) the values of mo, 1, jo, j, Ro, R, t, u, 1, (b) the principal terms in 
the expressions for T**, L*8, T«° (denoted by ‘x’). The terms not marked with ‘x’ are 
negligible in comparison those principal terms. It will be noted that there are two 
lines in the table for SS1, SS4, SS13, SS21, SFi, SF5, and three lines for SS19, 
SS20, SF4, SF8. This is because, in each case, k may have two or three values. 

For example, in the case of Type SS1, we have for T**, L8, 


Te = TP +T73?, Le = Li, 





while for T@°, 


To — T¢ (if k = 1), 
T®=T¢4+T7S (if k = 2). 
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THE AERODYNAMICS OF A RING AIRFOIL* 


BY 
H. J. STEWART 


California Institute of Technology 


Abstract. The downwash required to produce a given vorticity distribution is com- 
puted for a ring airfoil and the results are compared with the corresponding two- 
dimensional case. From this it appears that if the curvature of the chord plane is 
small, as is the case with normal amounts of dihedral, the effect of this curvature on 
the chordwise lift distribution of a wing is extremely small. If the radius of curvature 
is small compared to the chord, as it is near the vertex of a cranked wing, it is seen 
that this curvature may cause comparatively large changes in the lift distribution. 

1. Introduction. At the present time, the steady state two-dimensional airfoil the- 
ory is a highly developed subject; and, subject to the usual limitations of perfect 
fluid theory, solutions may be obtained with almost any desired degree of accuracy. 
The steady state three-dimensional airfoil theory is, however, in a much lower state 
of development. For most engineering problems the “lifting line” theory as developed 
by Prandtl and others is adequate to provide satisfactory results; however, for certain 
other problems, such as the flow near a wing tip, the effects of sweepback or of yaw, 
or the lift of a low aspect ratio wing, the lifting line theory cannot be used. At this 
time there have been only a fairly small number of solutions of finite wing problems 
in which lifting surfaces rather than lifting lines have been used and which may thus 
be used to throw light on these essentially more complicated problems. The best 
known of these lifting surface theories are those due to Blenk,! Kinner,? Krienes,* and 

Bollay.* As the number of such solutions is so 
limited almost any special solution involving a 

lifting surface is of interest. 
/ From an analytical viewpoint, probably the 
: Zz simplest lifting surface problem which has not 
\ \ yet been investigated is that of the axially sym- 
\ metric flow past a ring airfoil as shown in Fig. 1. 
\ \ This flow is especially simple as the vortex lines 
~~ in the lifting surface are circular rings and there 
| Lo are thus no trailing vortices. The particular 
+7 purpose of the present paper is to discuss the 
differences between this problem and the cor- 
Fic. 1. Ring wing in a uniform flow. responding two dimensional problem. In addi- 





. 





tion to its intrinsic interest in the theory of the 
“anti-drag” cowl, the ring airfoil problem possesses a general interest insofar as it 
demonstrates, at least qualitatively, some of the effects of dihedral on the lift dis- 
tribution of a wing. 





* Received Feb. 19, 1944. 
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2. The vector potential. The mathematical analysis of this problem may be con- 
veniently carried out by the method of the vector potential. Since the equation of 
continuity in an incompressible fluid is simply 


div q = 0, (1) 
the velocity vector g may be written as the curl of a vector potential A or 
q = curl A. (2) 


By the Helmholtz decomposition theorem the vector potential may be subjected to 
the restriction that 
div A = 0. (3) 

The differential equation for the determination of the vector potential is found by 
curling Eq. (2). This gives 

V7A = —curlg = —Q. (4) 
If the vorticity Q is a given function, this is a Poisson equation for the determination 
of the vector potential. The solution of this equation, which is well-known and may 
be obtained by the use of Green’s theorem, is 


1 dv 
Q— 


4r T1 


A (5) 


where the volume integral covers the entire region.where the vorticity exists and ri 
is the distance from the point at which the vorticity exists to the point P at which 
the vector potential is being computed. If the vorticity is in the form of a single 
vortex filament of strength [ then 
r 1 
A=— —ds, (6) 
4r 1) 
where ds is an infinitesimal distance vector along the vortex filament. If there are 
several vortex filaments the contribution from each one may be found by Eq. (6), and 
these results must then be summed to obtain the vector potential. 
3. The vector potential for a vortex ring. As the vortex filaments are all circles 
for the axially symmetric flow past a ring 





wing, the complete vector potential can easily Az 
be obtained if the vector potential of a single Ly p 
filament is known. For such a filament of de Ler 
strength [ and lying in the plane z=0 (see 
Fig. 2), it is obvious that the vector potential OZ 
ah see it “Nt oi 
is not a function of the meridian angle @, and — 
it may be calculated at points in the plane y 

+—-O 





6=0. Since ds is in the plane z=0, the vector 

potential can have no z-component. Further- ds 
more, by considering two vortex elements, 

one having the negative of the other’s @ co- 

ordinate, it is evident that the vector poten- 

tial can have no radial component. The vector potential has thus only the component 
A,» which is perpendicular to the meridian planes. By Eq. (6), this is 


Fic. 2. Vortex ring. 
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aI cos 6 dé 


Qe Jo [a? + 72 + 2? — 2ar cos 8]*/? 





Ag = 


With the vector potential expressed in this form as an elliptic integral, it is rather 
difficult to superimpose the vector potentials for a band of vorticity of radius a and 
of chord c in order to represent the ring wing. A much more convenient form can 
be obtained by the use of the Fourier integral. For an even function f(z), the Fourier 
integral theorem states that 


ae Lee 
(2) = — os kz: os kt dt>dk. 8 
f(z) — [cos 1 J 10 008 t (8) 


Since A, is an even function of z, it follows that 


al - 7 ad cos 6 dé 
Ay = — cos isl f cos kt if dt | dk. (9) 
rT? Jo 0 a [a? + r? — 2arcos@+ {2 ]1/2 


Since’ 








f cos (kt) dt = Ko(kz) (10) 
o [x?+ 2)1/2 sla iat 


the inner two integrals of Eq. (9) become, after inversion of the order of integration, 





I ={ Ko [kV a? + r? — 2ar cos 6| cos 6 dé. (11) 
0 


The addition theorem for the modified Bessel functions of the second kind (see Ref. 5, 
p. 74) states that 
To( ka) Ko(kr) + 3. I,( ka) K,(kr) cos n@ if r > a. 
Ko[kvV/a? + 7? — 2ar cos 6] = ya (12) 
To(kr)Ko(ka) + 2 >> I,(kr)K,(ka) cos nO if r < a. 


n=l 





Since the trigonometrical functions are orthogonal over the range 0S0S7, 
frli(ha) Ki (hr) if r>a. 


(13) 
lali(kr)Ki(ka) if r <a. 


The vector potential for the vortex ring in the outer range where r>a can thus be 


written as 


aI 7 
Ag = —f I,(ka)K,(kr) cos (kz)\dk = (r >a). (14) 
wr Jo 


For the inner range it is necessary to interchange the arguments of the two Bessel 
functions. 

4. The vector potential for a ring airfoil. A ring airfoil may be considered to be a 
system of ring vortices of radius a and distributed over the chord c from z= —c/2 
to s=c/2. If the strength of this vortex sheet is y(z)), then the vector potential 


forr>dais 


5 Grey, Mathews and MacRobert, Bessel functions, Macmillan and Co., London, 1931, p. 52. 
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a e/2 c) 

Ap = =f va) f 1,(ka)K,(kr) cos k(z — cae dzo. (15) 
wT J —¢/2 0 

If the vortex strength is known, Eq. (15), after inversion of the order of integration, 

can conveniently be used to compute the vector potential or the radial or axial velocity 

components, u, and u, respectively. From Eq. (2), 


OAg a. 
’ “«,>— aor (rAg). (16) 


Oz r or 





The radial velocity is of the most interest as it corresponds to the downwash velocity 
in the ordinary two-dimensional airfoil theory. The downwash at the ring airfoil, 


r=da, is 


i oo c/2 
uy = <f kI,( ka) K,(ka) if (zo) sin k(z — saldsoh dh (17) 


T c/2 


5. Comparison with two-dimensional flat plate airfoil. If the airfoil shape is given, 
the downwash u, is known, and Eq. (17) may be considered as an integral equation 
for the determination of the vortex strength y(z). It is, however, an integral equation 
of a difficult type. The importance of the curvature of the chord plane may be esti- 
mated by comparing the downwash for some given vortex distribution with the cor- 
responding two dimensional result. This process will be carried out for the vortex 


distribution 


G = 220 
(zo) = A s/ = ° (18) 

c + 220 
In the two-dimensional case, this vorticity distribution corresponds to a flat plate 
airfoil with its leading edge at 29 = —c/2. The downwash is then constant over the 


airfoil and equal to A/2. For this vorticity distribution it can easily be seen by use 
of the transformation 22) =c sin 8 that 


9 


fv sin k(z — 20)dzo = = Ac[olbke) sin kz + Ji(4kc) cos kz]. (19) 
The inaiinils velocity is thus 
u, = 4Aca f ” ply( ka) K (ka) [Jo($kc) sin kz + Ji(}kc) cos kz|dk. (20) 
0 
It is of interest to note that the two-dimensional result can be obtained directly 


from this by considering the limiting form as the radius of the ring becomes infinitely 


large; for 


lim { x7 1(x)Kx(x) } =}; (21) 
so the downwash in the two-dimensional case is given by 


uy = Ac f [albke sin (kz) + J(}hc) cos (kz) dk. (22) 
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TABLE [| 
Comparison of F(x) and F(x) 


F(x) F(x) | F(x) —F,(x) 





x 
0.1 0.0493 0.0132 0.0361 
0.5 0.2132 0.2000 0.0132 
1 0.3402 0.3637 —0.0235 
2 0.4450 0.4571 —0.0121 
3 0.4762 0.4800 —0.0038 
4 0.4873 0.4885 —0.0012 
5 0.4921 0.4926 —0.0005 





The first integral vanishes on the airfoil where 2? <c?/4 and the second is equal to 2/c 
on the airfoil (see Ref. 5, p. 65); so in the two-dimensional case, for this vorticity dis- 


tribution 
1 (2? S 2/4). (23) 


t|- 


ur = 


An exact evaluation of the integral of Eq. (20) is rather difficult; however, an ap- 
proximate evaluation, valid for large values of a/c, may be obtained quite easily. If 


F(x) = xI;(x) K(x), (24) 
a very close approximation to F(x) is given by 
x2/2 ; 
F(x) = ————__ (25) 
3/8 + x? 


It may be noted that the asymptotic expansions for F(x) and F,(x) are the same up 
through terms of order (x~?). It is shown in Table I and Fig. 3 that Fi(x) is a good 
approximation to F(x) even for small values of x. 






































0.6 
| | 
OF 
a2 
0 
0 / 2 x 3 4 5 


Since 
F\(x) -—3=- poten Ant (26) 


an approximate expression for Au, the difference between the ring airfoil downwash 
of Eq. (20) and the corresponding two-dimensional case is given by 
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Au=-— (3/32)4c f ol} ke) sin kz + J,(}kc) cos kz] _.. . (27) 
0 3/8 + a*k? 
Let \=kc/2, a=+/3/32c/a and B=2z2/c. Then 
Au = — sat f 00) sin BA + J:(A) cos Bd] =... (28) 
0 ? + a? 
On the airfoil where 8? <1, this gives (see Ref. 5, p. 78) 
Au = 3A [a cosh (a8)Ki(a) — a sinh (a8) Ko(a) — 1). (29) 


The ratio of the change in downwash to the two dimensional downwash is 
2(Au)/A. For a=0.02 and 0.20 corresponding to a/c=15.3 and 1.53 respectively, 
this ratio is given in Table 2 for the leading edge (@ = —1), the center of the airfoil 
(6 =0), and for the trailing edge (8 =1). 


TABLE 2 


Values of 2(Au)/A for the ring airfoil 








i 0.02 0.20 
a/c 15.3 1.53 
p= — 0.0009 si 0 .0451 
p=0 —0.0009 | —0.0448 
B= | 0.0023 |  —0.0963 





As the downwash velocity is determined by the slope of the camber line, the airfoil 
camber required to produce the lift distribution of Eq. (18) may be computed by in- 
tegrating the downwash velocity. The camber lines 
for a/c =1.53 and for the two dimensional case are 
shown for comparison in Fig. 4. 

6. Conclusions. From Table 2, it is apparent 
that the effects of the curvature of the chord plane 
of the ring airfoil are negligibly small if a/c =15.3 ee ee ee oe ee 
while they are fairly large for a/c=1.53. From vorticity. See Eq. (18). 

Fig. 4 it appears that the lift of a ring airfoil having 

a constant angle of attack across the chord would be somewhat more than that of 
the corresponding two dimensional airfoil and the lift is shifted away from the leading 
edge toward the center of the airfoil. 

It seems reasonable to suppose that the changes at any given section of the ring 
wing are caused primarily by the vortex elements near that section. These results 
may thus be applied in estimating the effects of the dihedral of a wing on the lift 
distribution over the wing’s surface. This indicates that if the curvature of the chord 
plane is small, as is normally the case, no appreciable changes in the lift distribution 
need be expected ; however, if the radius of curvature of the chord plane is of the same 
order as the chord, fairly large effects may be expected. This should be particularly 
noticeable near the vertex of a cranked wing. 
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THE MATHEMATICS OF WEIR FORMS* 


BY 
ALLEN P. COWGILL 
Syracuse University 


1. Introduction. This paper aims at making more readily available the results of 
a study of the mathematics of weir forms, a subject in Hydraulics to which higher 
mathematics can be applied. Section 2 covers the general application of Abel’s in- 
tegral equation to the forms of weirs by Brenke.! Section 3 deals with sectionally 
analytic weir forms, particularly the Stout-Sutro weir. The writer believes he has 
made the original application of Abel’s integral equation to this corrected weir form. 
Section 4 deals with cases when the quantity of flow can be expressed as a conver- 
gent series. 

2. Abel’s integral equation. One method of solution of the problem of weir forms 
involves Abel’s integral equation. The natural conditions found in the flow of water 
through weirs satisfy all the requirements of this integral equation, so it proves a 
superior mathematical tool in handling the general problem. In 1922 Brenke studied 
the problem of the weir form when the flow was proportional to some power of the 
depth. He made the original application of Abel’s integral equation. This equation 


has the form 





: ‘ (s)ds 
(2) = f eres, 0<r<1 (1) 
a (x cea s)* 
and its solution is, under certain conditions, 
sin Ar = ¢$'(s)ds 
*f(x) = f : (2) 
rg a (x—s) 


To obtain (2) from (1) use is made of two fundamental formulas, namely?* 


T 4 dx 
—_ =f —— 0<dv\<1 (3) 
sin Ar 2 (2 — x)**(x — s)* 


2 z U 1: Zz 1 zr , 1 
ff. th — ds = f f as nl (4) 
a w, (2 — x)'(x — s)* a (2 — x)*J, (x—s)* 


¢(s) is assumed to be continuous and have a continuous derivative in the closed 











* Received January 8, 1944. This paper constitutes part of a thesis submitted in partial fulfillment 
of the requirements for the degree of Master of Arts at the University of Nebraska. 
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1927, pp. 211, 229. 
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interval, a to b. Formula (4) is known as Dirichlet’s generalized formula.‘ Multiply (3) 
by #’(s)ds and integrate from a to 2, (a $z)), which gives 


T z z ’(s)d 
lo) —o@l= ff ee aii (5) 


sin Ar (z — x)!-*(x% — s)* 








If (4) is applied to the right hand member of (5), we have 


(2) — (a) = mes f € 1 "“#Oe , (6) 


z2—x)*J, (x-—s) 





Then, if é(a) =0 and if we replace the inner integral on the right of (6) by its value 
from (2), we see that (6) becomes (1). Hence (2) is a solution of (1). 

The weir is actually symmetrically constructed as in Fig. 3, but for purposes of 
the present calculation a half section is used (Fig. 1). Letting y=f(x) express the 


x 























Fie. 1. Fic. 2. 


distribution of width over depth, h the depth of flow, Ca the coefficient of discharge 
(approximately 0.6), and assuming that the quantity of flow is proportional to the 
mth power of the depth of stream, we have 


h 
caf [2¢(h — x) ]*/2f(x)dx = bh”, (7) 
0 
or, letting K =b/Ca(2g)'?, 


fo — x)'/?2f(x)dx = Kh™. 


Differentiating with respect to h, we have 


* f(x)dx 


———— = 2Kmh™. (8) 
0 (h — x)*/? 


This equation has the form of Abel’s integral equation. 
To find the equation of the weir form when the flow is bh, we have 


4W. A. Hurwitz, Note on certain iterated and multiple integrals, Annals of Math., 9, 183 (1907). 
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a sin a/2 (7 2Km(m — 1)h™? 
f(x) = lal f Bai: nts sd 
T 0 (x — h)}/? 
or 
; 2Km(m—1) 7* hm 
hin ee ay " 


7 o (x — h)'? 
By the use of Gamma Functions,5 


m—3/2 


f(2) = — wr, m 


IV 
to 


(10) 


Let be a positive integer 22. Then the Gamma Functions become simple products 


when m=n or m=n+3. When m=n, 


f(x) = — ———__— x32, n = 2. (11) 


When m=n-+H, 





’ 


1-3-5. +++ (2n+1) . 
~~ grt. eZ. (12) 
2"(n — 1)! 


3. Sectionally analytic weir forms. When m is equal to or greater than $ one 


gets continuous forms of weirs (Fig. 2). When m is greater than 3} and less than $ the 
weir forms have an infinite width at the bottom, the curve f(x) approaching the 
X-axis asymptotically. As this is impossible in practice, the necessary correction due 
to limiting the width of the weir furnishes an interesting mathematical problem which 
has been studied in the case where m =1. 





~ dh 
H 
| 

















ales Se 
{39 ; | 
A ———_—_—_—__—— WwW —— ees 
Fic. 3. Copy of Stout’s drawing in 1897. Fic. 4. 


The weir in which the flow is proportional to the depth is of engineering value. 
One of the first records of it is in an article by O. V. P. Stout. Approximate correction 
was made by circular openings (Fig. 3). A weir of this type was also constructed by 
Sutro and it is referred to in some texts as the Sutro weir. The modern way to correct 
5 F. S. Woods, Advanced calculus, Ginn & Co., 1926, p. 164. 

60. V. P. Stout, A new form of weir notch, Trans. of the Nebraska Engineering Society, 1, 13 (1897). 
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the Stout-Sutro weir is to start with a rectangular cross section of depth a and 
width w (Fig. 4). The upper section is then designed to give a flow proportional to 
the first power of the depth when the depth of flow equals or exceeds a. 

The calculations of E. A. Pratt’ by series solutions gave a mathematically correct 
form of weir where h2a. In this solution a rectangular section of depth a and width 
w is first assumed. Soundings are made with the zero point 3a from the.bottom, 


QO = bH = b(h + ja). 
The quantity of water discharged through the rectangular portion of the weir is 
Qo = $wK[(h + a)??? — h3/?], 
Therefore P 
QO = 4wK[(h + a)3/2 — f3/2] 4 2K f (h — x)'!*f(x)dx = b(h + §a). 
0 
As this equality must hold for h=0, 4wKa*/?=2ab and 6=2wKa"/?, so 


h 
f (h — x)"2f(x)dx = 2w[Zha'!? + a3/? — (h + a)*/? + hP/?], 


Instead of solving by the use of series, as Pratt did, one may differentiate with 
respect to h to put the equation in the form of Abel’s integral equation; thus 


h f(x)dx 
a _ 2w[a'/? — (h+a)'¥24 hil?], 
0 ;=— = = 


For the solution of Abel’s integral equation the right hand member must be a con- 
tinuous function, equal to zero when h=0. These conditions being satisfied, 











sin 4/2 = |— 3(h + a)? + §h-? 
y = f(x) = —— 2w f I dh 
, > 0 (x — kh)? 
w | z dh i dh ] 
= ' [ah au h?|}/2 0 [ax + h(x — a) — h?}/2 
Ww Ww ge=~-Z 
= —+ —sin™ , (13) 
2 T a+x 
or 
Qw x 1/2 
y= w- — tan-*(=) ‘ (14) 
vg a 
This solution can also be written 
a(w— y) 
x = a tan? —————_ - (15) 
2w 


In the design of the Stout-Sutro weir it is now necessary to choose an a for sub- 
stitution in the above formulas. One will generally know the average depth of flow 


7E. A. Pratt, Another proportional-flow weir, Sutro weir, Engr. News, 72, 462 (1914). 
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expected through the weir. It is felt to be better to keep the curve of (13) as close to 
the curve of the uncorrected weir derived from (10), y=2Kx~-/?/z, as possible. The 
scheme is to make the rectangular section of the Stout-Sutro weir have the same 
dimensions as if the uncorrected weir of (10) were corrected for the average depth of 
flow by the addition of a rectangular section at the bottom, below the Y-axis, to 
compensate for limiting its width to 2w. 

One substitutes y=w in the uncorrected formula (10) above and solves for x. 
This value and that of the 4 assumed to be average are substituted in 


ys 


“~ 


33/2 


h. , h — 2x in 2 ( ’ wh , 
mele poeeinand = a: ee 
> sin ; 3,in 1 + 1) r (hx x) 


v4 U 





(h — x)*/? = 0, 


which is solved for 7. One then makes a, the depth of the rectangular section, equal 
to x+r. It must be appreciated that (10) can be corrected for one depth of flow by 
the addition of a rectangular opening at the bottom, but would not be correct at 
any other depth of flow. Formula (13) is correct at any depth, H> a (Fig. 4). 

4. Series solutions of weir forms. We consider now the forms of weirs when the 
quantity of flow can be expressed as a convergent series in powers of h. Assume 
that the quantity of flow, Q(h), can be written 


n=2O 
Q(h) = >) ah", (16) 
n=0 
a convergent series not having a constant term, and assume the form of weir to be 
given by 


n=. 


f(x) = D> fal), (17) 


n=0 


each term of (17) giving rise to one term of (16). 
The general equation is 


h 
Ca(2g)* f (h — x)'2f(x)dx = Q(h). 
0 
Its solution will involve a series of integral equations of the form 
Ca(2g)2 f (h — x)*!*f,(x)dx = a,ht¢ n=0,1,2,---;a>}3 


which can be solved by the use of (10), giving 


Tin +a+t 1) 
fn(x) i Ca, oe ea my eer (18) 
T'(n + a — 3) 
where 


2 
C a(2gm)*/? 


Substitution of this in (17) gives the formal solution 


— Tin+t+at+i 
f(x) =C> >a, ee. gree, (19) 
0 I'(n+a-— }) 
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The first term of this series 


I'(a + 1) 
a TF 
I'(a — 3) 


fo(x) = Cao see . (20) 


will be discontinuous at x =0 if } <a<¥$ and continuous if a2. 

The series formed by all the terms after the first will converge and represent a 
continuous function of x. This may be proved as follows. By hypothesis the series 
>“""?a,h"t* converges since it is the series for Q(h), (16), with the first term 
omitted. Let 

Tin + a+ 1) 
= : ais 1 
64 ee ’ Qa 2 
T(n + a — 3) 
(n+ a)(n+a— 1)l(n+a-—1) 
= ; /  T(p+1) = pri), 
T'(n + a — 3) 


(n+ a)(n+a— I)c, 


where c,) =I'(n+a—1)/I'(n+a—}4) and 0<c, <1 since I'(p) increases monotoni- 
cally for p> 1.46. Now’ if the series }-"=a,h"* converges, so also will the series 











I 


N= N=00 n=00 
> chanh™**, Do netah™** and >> nc! anh™**. 
n=] n= 1 n=1 


n= 


also converges. In each case the function represented by the series is continuous. 
We have then the form of weir given by 


f(x) = fo(x) + g(x), (21) 
where fo(x) is given by (20) and 


But the series >-2=}°c,a,h"** is a simple combination of these three series, hence it 


n= 


g(x) = cy. CuBn a? te-8/2, 


n=1 


the quantities a, C and c, being as specified above. The solution f(x) is discontinuous 
at x=0 if 4<a<. It is continuous for a2 #. 


8 F. S. Woods, Advanced calculus, Ginn & Co., 1926, p. 47. 
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THE NUMERICAL SOLUTION OF LAPLACE’S AND 
POISSON’S EQUATIONS* 


BY 
DAVID MOSKOVITZ 


Carnegie Institute of Technology 


1. Introduction. A quite common method of solving numerically the Laplace dif- 
ferential equation 
avy #V 
mene ie enorme ee (1.1) 


Ox? =ay? 


with the boundary values of V prescribed on some contour I bounding a region R 
is to approximate V by the solution u of the Laplace difference equation: 


4u(x, y) = u(x+ h, y) + u(x — hh, y) + u(x, y+ bh) + u(x, y — hh). eee 


Briefly, the method of procedure, commonly called the Liebmann procedure,' is to 
cover the region R with a rectangular network of lines at distances h apart, and to as- 
sume values at the interior lattice points of this network. Using these assumed values 
and the known boundary values, we traverse the region R moving in some definite 
geometrical pattern from lattice point to lattice point, replacing the assumed values 
of u at each lattice point by the arithmetic average of the values of uw at the four 
neighboring lattice points. We then repeat the traverse moving in the same pattern 
to obtain a second improved value of u at each lattice point; and so on until a con- 
vergent stage is reached when the values of uw are no longer changed materially by 
continued traversing. 

The purpose of this paper is to present a process which yields precisely the con- 
vergent values of u obtained by infinitely many traverses of the region. In more 
precise language, if u; is the kth approximation of the value of wu after & traverses, 
our process yields the value u =lim;... ux. 

2. Notation and set up of the problem. Equation (1.2) can be transformed to 


4u(x, y) = u(x, y+ 1) + u(x, y — 1) + u(x + 1, y) + u(x — 1, y) (2) 


by a simple transformation, and we shall concern ourselves with the solution of equa- 
tion (2.1), with the values of u prescribed on the boundary lines 
x=0, x=n, y=0, and y=m 

of the rectangle R. 

Unless otherwise stated the numbers m and » are fixed positive integers, and 7 and 
j will be used as variable positive integers with the’range of values 

p21 2.+++,#=— I: g@i,2,---,m—I1. 

We shall denote the value of u at the point (7, 7) by u,;(z), and we desire to distinguish 


the known prescribed values of u on the boundary (which values are precisely the 





* Received Jan. 7, 1944. 
1 For a fuller description of the process and for references to the literature on the subject, see 


the paper by Shortley and Weller in J. App. Phys. 9, 334-348 (1938). 
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same as those of V on the boundary) from the unknown values of u at the interior 
lattice points. Accordingly, we denote the prescribed values of u on the boundaries 


as follows: 
by #,(0) at (0,7); by a(n) at (mn, 7); 
by wo(i) at (7,0); by tmn(i) at (i,m); 


and agree that 
u;(0) = u;(n) = u(t) = un(i) = 0 


wherever these terms appear in our equations. 
At each interior lattice point we can write ‘ 


4u;(i) = uj(i + 1) + ui — 1) + wil) + usd + o,(9), (2.2) 


where 
j(t) = 55,100 j(0) + 5:,n-10j(m) + 45,104) + 5;,m-1tm(9), (2.3) 
and 4,; (or 4;,;) is the Kronecker delta defined by 
i;=1, if i=; 6;=0, if t+). 
3. A special system of difference equations. We consider the solution of the sys- 
tem of equations 
Luj(i) = cuj(i) — ui + 1) — u(t — 1) = wgyi(t) + uja(d) + ¢,(9), (3.1) 
(¢=1,2,--+,n—1;7 =1,2,--+,m—1), 
in which c and the (m—1)(m—1) constants ¢,(¢) are prescribed.? We assume that 
u;(0) = uj;(n) = uo(i) = um(i) = 0, 


and seek the values of the (m—1)(m—1) unknowns 1;(1). 

The system (3.1) represents (m—1) difference equations in the (m—1) functions 
u(x), (7=1, 2, - - +, m—1), whose values are desired for integral values of the argu- 
ment x from 1 to (n—1). It can be readily shown that system (3.1) has a unique solu- 
tion, and of course this solution can be written down by Cramer’s Rule, but we shall 
give the solution in another form. . 

An immediate property of the operator L, defined by (3.1) is given by 


Le+at(t) - (L. + a)u(i) ’ (3.2) 


where a is any constant. We define the inverse operator Lz! and integral powers L? 
of the operator L, in the usual manner. 

From (3.2), we obtain 
[L. + a] = Life (3.3) 


and 
-1 


L.[Lere] oul Sone c+a (3.4) 
where 1 is used as the identity operator. The latter is established as follows: 
L.[Lij,] = [Leve — a] [Lope] = LeveL, — aL, = 1 — aLzyy. 


2 We assume that at least one of the ¢,(i) is different from zero; otherwise the system (3.1) has only 


the trivial solution u;(z) =0. 
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The solution of system (3.1) will be given symbolically in terms of the inverse 
operator L=', and the interpretation of the symbolic solution will depend on the con- 
stants D.(k) and X,,,(z) defined below. To apply the solution obtained to the solution 
of the Laplace difference system (2.2), we have merely to observe that (2.2) is a spe- 
cial case of (3.1) in which c=4 and ¢,(7) has the value given in (2.3). 

Let o, and p, be the roots of the characteristic equation 


c&é —-# —1=0. (3.5) 
We define D.(k) and X,,.(z) by 
ok _ pe 
uaa, (3. 6) 
Oc — Pe 
a 
he, g(t) = ——D.(n — i) when k Si, 
D.(n) 
(3.7) 
’ Dn — k) : 
he,k(4) = ————— D,(i) when k2i 
D.(n) 
The identities 
D.(n) = D.(i)D-(n — i+ 1) — Dit — 1)D-(n — 2) (3.8) 
and 
DAi + 1) = cD.(i) — D.(i — 1) (3.9) 


with D.(0) =0, D.(1) =1 are easily established. In terms of the operator Z., we may 
write (3.9) in the form 
L.D.(i) = 0. (3.10) 


Also if a is any constant, we have 


DsDas (7) = [Lore = a|De+a(i) = Le+eDe+a(4) aD. +.a(1); 


hence 
LeDesa(i) = — aDera(2). (3.11) 


We shall show that 
Lede, k(t) = Sit. (3.12) 


To establish (3.12) we have three cases to consider: 
1) When k Si—1, we have 
Lde,x(t) = ——L.D.(n — i) = 0, by (3.10). 


D.(n) 





2) When k21+1, we have 


D.(n — k) ; 
—— LD(i) =0, by (3.10). 
(2) 





Lede, x(t) = 
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3) When k =1, we have 
Lede,i(t) = Cre,s(t) — Ae, i(t + 1) — re i(t — 1) 
= [cD.(n — i)D.(i) — D-(n — i — 1)D,(i) — Di — 1)D.(n — i) ]/D-(n) 
= [D.(i){cD.(n — i) — D.(n — i — 1)} — D.(i — 1)D.(n — i) ]/D.(n) 
= [D.(i)D.(n — i+ 1) — D.(i — 1)D.(n — i)|/D.(n), by (3.9), 
= D.(n)/D-(n) = 1, by (3.8). 
Also useful is the relation 
Lede¢a,k(t) = 54,4 — Ode+a,k(4), (3.13) 
which is a consequence of (3.2) and (3.12). 
4. The special cases m =2 and m=3. As an introduction to the symbolic solution, 


we consider first the simplest case, m=2, in which case there is only one equation 
in the system (3.1), since by hypothesis uo(7) =u»(t) =0. This equation is 


L.u;(i) = $;(1). (4. 1) 
Its solution is given by 
n—1 
ui(i) = Dy re,n(i)pi (A). (4.2) 
kel 


We also write the solution of (4.1) in the symbolic form 
u(i) = Lz *[gx(i)], (4.3) 


where the expression on the right side of (4.3) is to be interpreted as being equal to 
the right side of (4.2). Explicitly, 


n—1 


L3*[$1(i) ] = Do re,e(i)bi (2). (4.4) 


To make actual use of the solution given in (4.2), it is necessary to have a table 
of values of the constants \,,.(4). In order not to complicate unduly the notation, 
the dependence of these constants on 7 has been omitted from our notation. A com- 
plete tabulation of these constants would require a great deal of space, since, with m 
fixed, they still depend on four parameters c, k, 1, and n. However, for the application 
to the solution of (2.2), we have c=4, and abridged usable tables requiring only 
tabulations for i=1, 2 and varying k and n are given in Tables 1 and 2 of §9. 

The entries in Table 1 give the multipliers to be applied to each of the boundary 
values in the calculation of ~;(1). As an example let us consider a 2 by 10 rectangle. 
Multiply each boundary value of the first column by the multiplier which appears 
opposite it in the column n=10; add these products, and the sum is the value of 
u;(1). Likewise, by interchanging the arguments (n—i) and 7, the same multipliers 
can be used to-calculate u:(9). Next using m:(1) and ~;(9) as known boundary values 
and the 2 by 8 rectangle which has the points (1, 1) and (1, 9) on its ends, use the 
multipliers in column n=8 to calculate u(2) and “3(8); and so on. 

The number of points at which the values of u are to be calculated by this process 
can be cut in half by using Table 2. With the entries in this table, using again a 2 by 
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10 rectangle, the values ~:(2) and (8) can be calculated using the multipliers in 
column =10; then in the 2 by 6 rectangle with points (2, 1) and (8, 1) as ends, 
and using multipliers in column n=6, calculate (4) and (6). The values 
u1(1), #:(3), - - + , u(7) can then be obtained by the Liebmann formula, each being 


the average of its four neighbors. For example, 
ur(S) = 4[ur(4) + u1(6) + @o(S) + a2(5)]. 
In the case m =3, system (3.1) reduces to the following: 
L.us(t) = uot) + gilt),  Letwa(t) = wilt) + ¢2(2), (4.5) 


with two unknown functions m and wu. 
Treating the operators in (4.5) as algebraic multipliers, and solving for u: and ue, 


we obtain symbolically 








m() = + 
(4.6) 
u2(i) = er he + E- 7 ont?) 


¢€ 
The following interpretation of the operators in the right members of (4.6) produces 
actual solutions. Write the operators of the right members in their partial fraction 
expansions, obtaining 


ace -4/ : + |= aust oh 
Cat "wt Les = —— 


1 -3/ 1 f- 1 ]|-3u2 —j- ]. 
a a?’ fae ii -_ 


c 














The last step in each line follows from (3.3). Consequently, the explicit solution of 
(4.5) is given by 


u(t) = . [an(i)bi(k) + Beli)pa(k)], a(t) = = [Bu(i)bs(k) + an(i)p2(h)]; (4.7) 


where 
oc;,(4) = 4 [r--1,4(4) + e+1,2(2) |, Bx(4) = 4 [re-1,2(4) - Ae41, 4 (2) |. (4.8) 


To verify that (4.7) is the solution of (4.5), we use 
Leax(t) = dix + Bx(4), L.Bx(t) = ax(i), 
which follow from (4.8) and (3.13). Consequently, 


n—1 


Lui) = >> { [Lecex(i) ]6:(k) + [L.Bx(2) ]o2(&) } 


kal 


n—1 


= F buds) + ¥ [Bldox(k) + an(ib2(8)] 


k=1 


= $:(t) + u2(%). 
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Tables 3 and 4 of §9 are useful in calculating values of u;(1) and (2) for the case 
m =3, again using c=4 which is appropriate for the Laplace equation. 

These tables are used in a similar manner as Tables 1 and 2. However, for each 
value of , there are two numbers side by side, and also two boundary values in the 
first column. The interpretation is to multiply the left one of the two boundary values 
of each line by the left one of the two multipliers of the same line in the appropriate 
n column. 

For a given , calculate u;(1) using multipliers of Table 3; then latorchenaiog the 
subscripts 1 and 2 and the subscripts 0 and 3, calculate u2(2) using multipliers of 
Table 4. The value (1) is then obtained by the Liebmann formula 


u2(1) = 3[2(0) + as(1) + ui(1) + u2(2)]. 


As in the case m=2, we can also calculate from the end x = of the rectangle. 

5. The general case. We now consider the general system (3.1) for any value 
of m, and show that its solution can be written symbolically in terms of the poly- 
nomial operators defined by 


Py, = 0, P, = I, P, = L.Pi-i1 — Pi-2 for >t (5.1) 


The operator P; is a polynomial in L, of the (k—1)st degree, which is precisely the 
same function of L. as D.(k) is of c. From this observation, the following analogue of 
(3.8) can be seen to be valid, 


Pa - PP. = Py1P 2+ (5.2) 


By solving the system (3.1) treating L, as an algebraic multiplier, we obtain sym- 
bolically 


i Pp P;Pm-t 
uj(i) = >> ——— i + ee 


kel ?. kejt+t «=P m 


——— 6), (U=1,2,---,8—2), (5.3) 


in which the operators of the right member are to be interpreted similarly to those 
of (4.6); that is they are to be expanded into partial fractions. Since each of the opera- 
tor coefficients in the right member of (5.3) is a proper fraction, their expansions will 


have the form 


PP he 
IRs = ES. Fas ore. (5.4) 
Fu lel | a lel 


where 41, d2, - - * , @m—, are the roots of P,,(£) =0 and the numbers b,(j, &) are uniquely 
determined. The actual solution of (3.1) is thus given in terms of the operators Lz", 
which have the meaning given in (4.4). 

That the foregoing actually yields the solution of (3.1) can be established by oper- 
ating on each side of (5.3) with the operator L, and using the relations (5.1) and (5.2). 

6. Application to solutions of Laplace’s equation. We have already observed that 
(2.2) is a special case of the system (3.1) in which c=4 and ¢,(7) has the value given 
in (2.3). To apply the preceding results, it becomes necessary to calculate the coefh- 
cients of ¢;(#) in the solution of (3.1) for various values of m. The polynomial opera- 
tors P,, must be factored, and the operators appearing in (5.3) must be written in the 
more useful form (5.4). 
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We have already considered in detail the cases m=2, m=3 in §4. We now par- 
ticularize the general solution of §5 to the case m =4. 
From (5.3) we have symbolically: 


















































P;P3 ii P,P» 4 P,P; ‘ 
1 = “ae 2 ’ 
41 P, 1 P, co) P, 3 
i OO ak SP el (6.1) 
Uy = 2 — A ). 
ee _ 
P,P, * P,P» + P,P; 
U3; = P, o1 P, 2 P, ak 
where 
P,P; 1 if 1 1 3 ; 
= = —————— = $/ ——__- + - = — — | = BL +Lapva— Lz" 
; 8-8, “hast £4.46. & i hin2te",| 
P,P, i. 1 1 1 - 4 
i a! on Oe 
P, iS «= 22. 2/2LL.—~vV2 Le+ V2 2/2 | (6.2) 
(6.2 
ne | ee =| t [Let ve +L vet 2L2"], a 
——T = — ~- — = c+ v2 
~ £2, “tht £400. ks) * | 
_ . | _— | gL vet Lave] 
; fae Liew £44); °° ° 
The explicit form of the solution is given by 
uy(i) = z [Qe(i)o1(k) + Su()bo(k) + Ri(ios(&)], 
k=1 
n—1 
u2(i) = D> [Se(i) {1(k) + $3(k)} + Te(a)ba(A)], (6.3) 
k=] 
n—1 
us(i) = >> [Re(i)oi(k) + Sx(i)b2(k) + Ox(i)63(2) ], 
k=1 
where 
Qx(i) = 2 [dc—va,0(i) + Aer v3,0(4) + 2Xe,e(4) J) ) 
Ri(i) = 4 [he—va,n(i) + Ae v3,4(4) — 2e,x(i) J, 
Eee : (6.4) 
Si(i) = mW, e~V3,k(4) — dey v3,4(4) ], 
T(t) = 4 [d-—v3,4(4) + e+ V3, 4(4) |, 
and 


1(k) = 5x,1%1(0) + bx,n-101(m) + tdo(R), 
go(k) = 5¢,1%2(0) + 54,n-1%2(m), (6.5) 
b3(k) = Sx,1%3(0) + Sk,n-143(m) + a(R). | 


Multipliers for calculating u2(1), u2(2), and u:(2) are given in Tables 5, 6, and 7. 
For a 4 by ” rectangle, calculate u2(1) using Table 5, u:(2) using Table 7, w3(2) using 
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Table 7 with subscripts 0 and 4 and subscripts 1 and 3 interchanged. Then calculate 
ui(1), us(1) by the Liebmann Formula. Then using known values u(1), u2(1), u3(1) 
and treating them as known boundary values for the rectangle one unit shorter, calcu- 
late u2(3) using Table 6; u2(2) can then be obtained by use of the Liebmann Formula. 
Then calculate :(4), us3(4) using Table 7; and so on. Of course, as in the case m =2 
and m=3, we can also calculate from the end x =m of the rectangle. By this process 
the value at only every other point on each line is calculated by the use of the tabu- 
lar values. . 

For values of m larger than 4, it is only convenient to use composite values of m. 
It can be easily shown that when m is composite having gq for one factor then P,, con- 
tains P, as a factor. However, for larger values of m, the necessary tables occupy more 
space and require a tremendous amount of time in their preparation. Theoretically the 
complete solution for any value of m is given by (5.3), but practically it is sufficient to 
use tables with m=4. A given rectangle can be covered by a lattice 4 units wide and 
the values of the function u at the interior lattice points calculated by the methods 
indicated. When the values of u at these lattice points are obtained, each rectangle can 
then be subdivided again by a finer network and the first calculated values are then 
used as approximate values for the finer network, and can in turn be improved either 
by traversing or by the use of the tables given here. 

7. The Poisson Equation. The preceding methods and results can be extended 
with slight modification to apply to the numerical solution of the more general Pois- 
son equation 

aV @V 

an tap  h™ ¥), (7.1) 
in which F(x, y) is defined in the interior of the region R. The approximating differ- 
ence equation in this case is 


4u(x, y) = u(x +h, y) + u(x — h, y) + u(x, y + A) + u(x, y — h) — WPF (x, y). (7.2) 
Employing our former notation, we can write at each interior lattice point 
4uj(i) = uit 1) + u(t — 1) + ugyalad) + uja(t) — WF (i) + 8:10 ;(0) + 5,10 ;(n) 

+ 5 ;,1%o( i) + 5 j,m—1tm(1), (i - 1, 2, PS gan te - t, 2, ‘S: + 1), (7.3) 


and again consider 
uj(0) = u;(n) = uo(t) = un(i) = 0. 


The system (7.3) is again a special case of system (3.1) in which the known func- 
tion ¢,(7) is now given by 
(i) = 55,10 ;(0) + 5;n-100j(m) + 5;,1000(7) + 5j,m—1m(t) — A°F ;(i). (7.4) 
Consequently, the general solution given in §5 applies at once. 
In the case m =2, for example, we have 


ui(i) = s e,e(4) [5x,19%1(0) + Se,n—1t41(m) + ao(k) + ti2(k) — hF (k) |. (7.5) 


To apply Tables 1 and 2 to obtain the values of (1) and m(2), first calculate 
h*F\(k), k=1, 2, - ++ ,n—1, at each interior lattice point. Then apply to these values 
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the same multipliers as are applied to #(k) and #(k). If Table 2 is used to calculate 
u;(2), the value u:(1) can be obtained from the associated Liebmann equation 


uy(1) = 3[m(2) + 22,(0) + ao(1) + (1) — WF ,(1)]. (7.6) 


For the case m =3, multiply h?Fi(k) and h?F2(k), respectively, by the same multi- 


pliers as are used for #o(k) and @;3(k). 
For the case m =4, additional tables are required. The explicit solution in this case 


is formally the same as that given in (6.3), but ¢:(k), ¢2(k), and ¢3(k) have the fol- 
lowing values: 


$1(R) ” 5;.,1%1(0) + 5% n—1041(m) + tlo( k) ea WF \(k), 
o2(k) = 5x,14%42(0) + 5k,n—1M2(n) mae h’F.(k), (7.7) 
3(k) = 5x,1%3(0) + 5¢,n-1%3(m) + 4(k) — h?F 3(k). 


The multipliers Q,(2) and R,(2) appear in Table 7 and are respectively those 
multipliers applied to #o(k) and @(k) in the calculation of u;(2). The multipliers S,(1) 
appear in Table 5 and are those multipliers applied to both #o(k) and #(k) in the 
calculation of u2(1). The multipliers S,(2) are those multipliers in Table 6 which are 
applied to both #(k) and #,(k) in the calculation of u2(2). 

The multipliers 7,(1) and 7,(2) which must be applied to h?F,(k) in the calcula- 
tion of u2(1) and “2(2) appear in Tables 8 and 9 respectively. 

8. Irregular cases and non-rectangular boundaries. The preceding solutions both 
in the case of the Laplace equation and the Poisson equation apply only to rectangu- 
lar boundaries whose dimensions are integral multiples of the lattice unit h. To apply 
Tables 5, 6, and 7 in the solution of the Laplace equation for a rectangular boundary, 
divide the smaller dimension of the rectangle by 4 to obtain the lattice unit h. If the 
longer dimension of the rectangle is an integral multiple of 4, the process here outlined 
for the solution applies directly. We call this the regular case. If the longer dimension 
is not an integral multiple of hk, we call this the irregular case, and a modification of 
the process here outlined is required. We need an analogue of the Liebmann Formula 
to express the value of a harmonic function approximately in terms of the values of 
the function at four non-equidistant neighbors. 

Let H(x, y) be an arbitrary harmonic function whose value Ho at (xo, yo) is to be 
expressed approximately in terms of Hi, He, H3, Hy which are the values of H(x, y) 
at the points (xo+rih, yo), (xo—reh, yo), (xo, Yotrsh), (xo, Yo—rat) where 11, 72, 73, 14, 
and f/ are positive. When 1 =72=r3=7,=1, and H(x, y) is approximated by its T. S. 
(Taylor Series) expansion about (xo, yo) up to terms including those of the third de- 
gree in h, Hp is found to satisfy the Liebmann equation 


Hy = 3(Mi1 + H2+ Hs + A). (8.1) 


When 11, 72, 73, and 74 are not equal, and H(x, y) is approximated by its T.S. up to 
terms including those of the second degree in h, we find? 
Ho = a;Hi + 2H, + asHs + aH (8.2) 


where 


3 This relation is also given by Shortley and Weller, ibid.; but their method of derivation is slightly 


different from ours. 
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a, = rorsta/ (ti + r2)(rite + ara), = a = Har ata/ (11 + 12) (rive + 774), 
8.3 
@3 = ryrora/ (13 + 14) (rire + 1374), O, = ryrors/ (13 + 14) (rite + 374). 
If H(x, y) is approximated by its T.S. up to terms of the first degree in h, we find 
4 
Ho = > b.H, with by = 1/7! + rz? 4+ 73? 4+ 17). (8.4) 


k=1 


Either (8.2) or (8.4) expresses Hy as a weighted average of its four neighbors, and 
although (8.4) is easier to use, presumably (8.2) gives a better approximation to the 
value of Ho. However in the application to the irregular case of the rectangle, we 
require only the simpler forms to which (8.3) and (8.4) reduce when three of the 
r. (k=1, 2, 3, 4) are equal to unity. 

To determine the values of a harmonic function at the interior lattice points of a 
rectangle whose dimensions are 4h by (+1r)h where n is an integer and 0<r<1, 
let x, y, 2, etc. be the values of the harmonic function at the points indicated in Fig. 1. 
By means of Tables 5, 6, 7 we can express u, v, and w as linear functions of x, y, 2 
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and the boundary values. Then by means of (8.2) or (8.4) or any other approximation 
method determine x, y, and z in terms of u, v, w, and the boundary values B, 
(k =1, 2, 3, 4, 5) at the indicated points of the figure. This process leads to three linear 
algebraic equations for the determination of x, y, and z. When these values are de- 
termined, the values of the harmonic function at the other interior lattice points can 
be obtained as the problem is now reduced to the regular case. 

A similar method can of course be used for the Poisson equation. In this case the 
analogue of (8.2) is 


Ho = aH, + a2H2 + a3H3 + aH, — aolo (8.5) 
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where 
F P 1 Vol 374 
a = +h? ——-—— }, 
rire + 1304 (8.6) 


Fy denotes F(xo, yo), and a1, @2, @3, ds are given by (8.3). 

The foregoing method applies equally well if the top and bottom boundaries of 
the figure are not straight. We postpone for a later paper the procedure which can 
be applied for non-rectangular boundaries in general. In this subsequent paper we 
shall also show the application of our methods to an extension of the Liebmann 
process in which the values at certain of the interior lattice points are calculated 
from the arithmetic average of the values at their four normal neighbors while the val- 
ues at the other lattice points are calculated from the values at their four diagonal 
neighbors. 

9. Tables. The entries in the following tables were rounded off to four decimal 
places from calculations carried out to a higher number of decimal places. In the 
tables, the decimal points are not printed but are to be understood to be present just 
before the first digit. In the compilation of these tables, the values of D.(k), defined 
in (3.6), were required for the values c=4, 3, 5, 4—1/2, and 4+-4/2. These were cal- 
culated from the recurrence relation (3.9), and are integers only when c is an integer. 

The values \4,,(1) and d4,,(2), defined in (3.7) are the entries in Tables 1 and 2, 
respectively. The entries in Tables 3 and 4 are the values a;(z) and 8;(7) fori =1 and 2, 
respectively; these were calculated from \3,;(z) and \5,x.(4) for i=1 and 2 using the 
relations (4.8). The entries of Tables 5 to 9 were calculated from the relations (6.4). 

As one convenient check on the accuracy of Tables 1 to 7, the sum of the multi- 
pliers used on all of the boundary values must be equal to unity; this check was 
applied. The author would be happy to know that no errors were made in the many 
calculations required in the preparation of these tables. 
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Boundary ms 
waliion n=3 n=4 2n=5 | n=6 n=7 n=8 n=9 
(0) | 2667 2679 | 2679 | 2679 | 2679 2679 2679 

tio( 1) +12(1) 2667 2679 | 2679 2679 2679 2679 2679 

#io(2) +142(2) 0667 0714 0718 0718 0718 0718 0718 

tio(3) + 122(3) 0179 0191 0192 0192 0192 0192 

tio(4) +a2(4) 0048 0051 0052 0052 0052 

fio(5) +142(5) 0013 0014 0014 0014 

tio(6) + t2(6) 0003 0004 0004 

io(7) +12(7) 0001 0001 

tio(8) + 12(8) 0000 
ti,(n) 0667 0179 | 0048 | 0013 | 0003 0001 | 0000 

TABLE 2.—To calculate u;(2); m=2 

Boundary 
: is n=4 n=5 n=6 n=] n=8 n=9 n=10 
values 
,(0) 0714 0718 | 0718 | 0718 0718 0718 | 0718 

tio(1) +22(1) 0714 0718 0718 0718 0718 0718 0718 

fio(2) +1%2(2) 2857 2871 2872 2872 2872 2872 2872 

io(3)+0,(3) | 0714 0766 0769 0770 0770 0770 0770 

tio(4)+m(4) | 0191 0205 0206 0206 0206 0206 

fio(5)+0,(5) | 0051 0055 0055 0055 0055 

ti(6)+12(6) | 0014 0015 0015 0015 

fio(7)+a2(7) | 0004 0004 0004 

fio(8)+12(8) | 0001 0001 

tio(9) +-222(9) 0000 
ti,(n) | 0714 o191 | 0051 | 0014 | 0004 0001 | 0000 

TABLE 3.—To calculated u(1); m=3 

B lary 
pr iti n=2 n=3 n=4 n=§5 n=6 
values 

2,(0) | a(0) | 2667 | 0667 | 2917 | 0833 | 2948 | 0861 | 2953 | 0866 | 2953 | 0866 
o(1) | (1) | 2667 2917 | 0833 | 2948 | 0861 | 2953 | 0866 | 2953 | 0866 
io(2) | (2) | 0833 | 0417 | 0932 | 0497 | 0945 | 0509 | 0947 | 0511 
o(3) | (3) | 0282 | 0195 | 0318 | 0227 | 0323 | 0232 
o(4) | a(4) | 0100 | 0082 | 0114 | 0095 
tio(5) | 15(5) 0037 | 0033 

tin(n) | 2667 | 0667 | 0833 | 0417 | 0282 | 0195 | 0100 | 0082 | 0037 | 0033 


i, (n) 
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TABLE 3.—(continued) 

Boundary | = | > 
values <spigd n=8 | n=9 n=11 
(0) | m(0) | 2953 | 0866 | 2953 | 0866 | 2953 | 0866 | 2953 - 0866 | 2953 | 0866 
ao(1) | as(1) | 2953 | 0866 | 2953 | 0866 | 2953 | 0866 | {2953 | 0866 | 2953 | 0866 
fio(2) | #5(2) | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 |40947 | 0512 | 0947 | 0512 
fio(3) | (3) | 0324 | 0233 | 0324 | 0233 | 0324 | 0233 [$0324 | 0233 | 0324 | 0233 
tio(4) | a9(4) | 0116 | 0097 | 0116 | 0097 | 0116 | 0097 |f0116 | 0097 | 0116 | 0097 
tio(5) | a(S) | 0042 | 0038 | 0043 | 0039 | 0043 | 0039 | 0043 | 0039 | 0043 | 0039 
(6) | (6) | 0014 | 0013 | 0016 | 0015 | 0016 | 0015 | 0016 | 0015 | O016 | o015 
tio(7) | as(7) | 0005 | 0005 | 0006 | 0006 | 0006 | 0006 | 0006 | 0006 
fi(8) | a3(8) | 0002 | 0002 | 0002 | 0002 | 0002 | 0002 
(9) | t5(9) | | | | | 0001 | 0001 | 0001 | O001 
fio(10) a4(10) | | | | | | 0000 | 0000 
ti(n) | a(n) | 0014 | 0013 | 0005 | 0002 | 0002 | 0001 | 0001 | 0000 | 0000 
TABLE 4.—To calc ulate whit m=3 

Boundary 
oniaee | n=3 n=4 | n=5 | n=6 n=7 
(0) | (0) » | 0833 | 0417 | 0932 | 0497 | 0945 | 0509 | 0947 | 0511 | 0947 | osi2 
io(1) | as(1) | 0833 | 0417 | 0932 | 0497 | 0945 | 0509 | 0947 | 0511 | 0947 | 0512 
fio(2) | as(2) | 2916 | 0834 | 3230 | 1056 | 3271 | 1093 | 3277 | 1098 | 3277 | 1099 
o(3) | (3) | | 0932 | 0497 | 1045 | 0591 | 1061 | 0606 | 1063 | 0608 
fio(4) | as(4) | | | 0318 | 0227 | 0360 | 0265 | 0366 | 0271 
io(5) | 19(5) 0114 | 0095 | 0129 | 0109 
1o(6) | 2(6) | | | 0042 | 0038 
a(n) | a(n) | 2916 | 0834 | 0932 | o497 | 0318 | 0227 | 0114 | 0095 | 0042 | 0038 
TABLE 4. —(continued) 
Boundary | 
alien n=8 | n=9 | n=10 n=11 n2=12 
(0) | @(0) | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 
fio(1) | (1) | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 | 0947 | 0512 
o(2) | aa(2) | 3277 | 1099 | 3277 | 1099 | 3277 | 1099 | 3277 | 1099 | 3277 | 1099 
Ho(3) | a(3) | 1063 | 0609 | 1063 | 0609 | 1063 | 0609 | 1063 | 0609 | 1063 | 0609 
ao(4) | a(4) | 0367 | 0272 | 0367 | 0272 | 0367 | 0272 | 0367 | 0272 | 0367 | 0272 
(5) | a(5) | 0131 | 0112 | 0132 | 0112 | 0132 | O112 | 0132 | 0112 | 0132 | 0112 
fio(6) | a(6) | 0048 | 0044 | 0049 | 0044 | 0049 | 0045 | 0049 | 0045 | 0049 | 0045 
(7) | #(7) | 0016 | 0015 | 0018 | 0017 | 0018 | 0017 | 0018 | 0017 | 0018 | 0017 
1o(8) | (8) | 0006 | 0006 | 0007 | 0007 | 0007 | 0007 | 0007 | 0007 
(9) | (9) | | | 0002 | 0002 | 0003 | 0003 | 0003 | 0003 
fio(10) | (10) | | 0001 | 0001 | 0001 | 0001 
fio(11) | #(11) | | 0000 | 0000 
a(n) | a(n) | 0016 | 0015 | 0006 | 0006 | 0002 | 0002 | 0001 | 0001 | 0000 | 0000 
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TABLE 5.—To calculate u2(1); m=4 











Boundary 


aliens n=3 | n=4 | n=5 | n=6 | n=7 | n=8 | n=9 














ti2(0) 2857 | 3230 | 3304 | 3320 | 3323 | 3324 | 3324 | 3324 | 3324 | 3324 | 3324 | 3324 





4, (0) +243(0) | 0714 | 0932 0982 | 0993 | 0996 | 0997 | 0997 | 0997 | 0997 | 0997 | 0997 | 0997 





| 
| 
tho(1) +-144(1) ‘as 0932 | 0982 | 0993 | 0996 | 0997 | 0997 | 0997 | 0997 | 0997 | 0997 | 0997 





fio(2) +14,(2) 0497 | 0625 | 0654 | 0661 | 0662 | 0663 | 0663 | 0663 | 0663 | 0663 | 0663 
io(3) +2,(3) | 0268 | 0332 | 0346 | 0349 | 0350 | 0350 | 0350 | 0350 | 0350 | 0350 
fio(4) +24(4) | 0133 | 0164 | 0171 | 0172 | 0173 | 0173 | 0173 | 0173 | 0173 
fio(5) +004(5) | 0064 | 0079 | 0082 | 0083 | 0083 | 0083 | 0083 | 0083 
io(6) +14(6) | 0031 | 0038 | 0039 | 0040 | 0040 | 0040 | 0040 
tio(7) +244(7) | 0015 | 0018 | 0019 | 0019 | 0019 | 0019 
tio(8) +244(8) | 0007 | 0008 | 0009 | 0009 | 0009 
(9) +2,(9) | 0003 | 0004 | 0004 | 0004 
to(10) +-1,(10)| 0002 | 0002 | 0002 
tho(11) +2%,(11)} 0001 | 0001 
fo(12) +124(12)} 0000 











ti, (n) +43(n) | o714 | 0497 0268 | 0133 | 0064 | 0031 | 0015 | 0007 0003 | 0002 0001 | 0000 



































a(n) | 2857 | 1056 | 0446 | 0226 | 0093 | 0044 | 0021 0010 | 000s | 0002 0001 | 0000 





TABLE 6.—To calculate u2(2); m=4 

















Boundary 
po n= 3] nad n=5in=6|n=7|n=8|n=9 oe Bee ee See ee oes 
#,(0) | 1086 1250| 1292) 1301] 1303] 1304| 1304| 1304 1304) 1304] 1304] 1304! 1304 
4,(0) +2,(0) | 0497) 0625| 0654| 0661! 0662| 0663| 0663| 0663| 0663| 0663| 0663! 0663| 0663 





tho(1) +-22,(1) 
tho(2) +1%,(2) 
tho(3) +-1%4(3) 0625) 0788] 0825) 0833) 0835) 0835) 0835) 0835) 0835) 0835) 0835) 0835 


| 0497 0625| 0654! 0661] 0662| 0663| 0663] 0663) 0663| 0663] 0663| 0663| 0663 
tho(4) + 144(4) | 0332) 0410) 0428) 0432) 0433| 0433] 0433) 0433| 0433] 0433] 0433 
| 


0932) 1250) 1325) 1342) 1346) 1347| 1347) 1347) 1347) 1347) 1347) 1347) 1347 


0164} 0202) 0210} 0212) 0212) 0212) 0212; 0212) 0212) 0212 

















tio(5) +444(5) 

tio(6) + 124(6) 0079} 0097) 0101; 0102) 0102} 0102) 0102) 0102; 0102 
tho(7) +1%4(7) 0038; 0046; 0048) 0048) 0049) 0049) 0049) 0049 
tio(8) +124(8) 0018} 0022! 0023) 0023) 0023) 0023) 0023 
Ho(9) + 224(9) 0008} 0010; 0011) 0011) 0011) 0011 
tio( 10) +24,(10) 0004} 0005} 0005) 0005) 0005 
tio( 11) +-22,(11) 0002) 0002; 0002; 0002 
tio(12) +-2%,(12) 0001) 0001) 0001 
#io( 13) +-22,(13) 0000} 0001 
to( 14) +-2%,(14) 0000 
2, (2) +22;(2) 0932) 0625) 0332) 0164) 0079) 0038) 0018 anne 0004; 0002; 0001) 0000) 0000 











a(n) | 3230 1250 0539] 0245) 0114 0054) 0025 0012, 0006 0003| 0001; 0001; 0000 
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TABLE 7.—To calculate u;(2); m=4 

Boundary = = — | _ a * 

wales n=; n=4 | n= 5 | n=6 n=7 } n=8 

(0) 0497 0625 | 0654 | _ 0661 0662 =| 0663 
(0) | as(0) | 0861 | 0195 | 0982 | 0268 | 1005 | 0287 | 1010 | 0292 | 1011 | 0293 | 1011 | 0293 
o(1) | a%(1) | 0861 | 0195 | 0982 | 0268 | 1005 | 0287 | 1010 | 0292 | 1011 | 0293 1011 | 0293 
fio(2) | (2) | 2948 | 0282 | 3304 | 0446 | 3365 | 0494 3377 | 0506 | 3380 | 0508 | 3381 | 0509 
fio(3) | 24(3) | 0982 | 0268 | 1129 | 0364 | 1158 | 0389 | 1164 | 0394 | 1165 | 0396 
fio(4) | 14(4) | | 0365 | 0174 | 0429 | 0224 | 0442 | 0236 | 0445 | 0239 
to(5) | 24(5) | 0148 | 0097 | 0177 | 0122 | 0183 | 0128 
tio(6) | a4(6) | | | 0064 | 0050 | 0077 | 0062 
fio(7) | ts(7) | 0029 | 0025 

! 

,(n) | ais(n) | 2948 | 0282 | 0982 | 0268 | 0365 | 0174 | 0148 | 0097 | 0064 | 0050 | 0029 | 0025 
tio(n) 0932 0625 0332 0164 0079 0038 
TABLE 7.—(continued) 

Boundary | | 

‘silane n=9 n=10 n=11 n=12 n=13 n=14 

ti,(0) 0663 0663 0663 0663 «=| = 0663 0663 
said | _——|—__—— ier : a ae a 
(0) | #(0) | 1011 | 0293 | 1011 | 0293 | 1011 | 0293 | 1011 | 0293 | 1011 0293 | 1011 | 0293 
ao(1) | a(1) | 1011 | 0293 | 1011 | 0293 | 1011 | 0293 | 1011 | 0293 | 1011 | 0293 | 1011 | 0293 
tio(2) | i%4(2) | 3381 | 0509 | 3381 | 0509 | 3381 | 0509 | 3381 | 0509 | 3381 | 0509 | 3381 | 0509 
fio(3) | (3) | 1166 | 0396 1166 | 0396 | 1166 | 0396 | 1166 | 0396 | 1166 | 0396 | 1166 | 0396 
tio(4) | a4(4) | 0446 | 0240 | 0446 | 0240 | 0446 | 0240 | 0446 | 0240 | 0446 | 0240 | 0446 | 0240 
tio(5) | (5) | 0184 | 0129 | 0185 | 0129 | 0185 | 0130 | 0185 | 0130 | 0185 | 0130 | 0185 | 0130 
io(6) | %4(6) | 0080 | 0065 | 0081 | 0066 | 0081 | 0066 | 0081 | 0066 | 0081 | 0066 | 0081 | 0066 
tio(7) | (7) | 0035 | 0031 | 0036 | 0032 | 0036 | 0032 | 0037 | 0033 | 0037 | 0033 | 0037 | 0033 
fio(8) | %4(8) | 0013 | 0012 | 0016 | 0015 | 0017 | 0016 | 0017 | 0016 | 0017 | 0016 | 0017 | 0016 
tio(9) | (9) | 0006 | 0006 | 0007 | 0007 | 0008 | 0008 | 0008 | 0008 | 0008 | 0008 
fi9(10)| 1,(10) | | 0003 | 0003 | 0003 | 0003 | 0004 | 0004 | 0004 | 0004 
fio(11)| 124(11)! | 0001 | 0001 | 0002 | 0002 | 0002 | 0002 
fio(12)} 24(12) 0001 | 0001 | 0001 | 0001 
fio(13)) 24(13) | 0000 | 0000 
h(n) | a(n) | 0013 | 0012 | 0006 | 0006 | 0003 | 0003 | 0001 | 0001 | 0001 | 0001 | 0000 | 0000 

0018 0008 0004 0002 | 0001 0000 


ii,(n) 
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TABLE 8.—Values of 7;(1) 





7 Se me 






































| | | | | | 0001 | 0001 
14 | | | | | | | 0000 


k | n=2 | n=3 n=4|n=5|n=6|n=7|n=8 | n=9 |n=10|\n=11/n=12|n=13 n2=14 
1 | 2857 | 3230 | 3304 | 3320 | 3323 | 3324 | 3324 | 3324 | 3324 | 3324 | 3324 | 3324 | 3324 
2 1056 | 1250 | 1292 | 1301 | 1303 | 1304 1304 | 1304 | 1304 | 1304 | 1304 | 1304 
3 0446 | 0539 | 0560 | 0564 | 0565 | 0565 | 0566 | 0566 | 0566 | 0566 | 0566 
4 0201 | 0245 | 0255 | 0257 | 0258 | 0258 | 0258 | 0258 | 0258 | 0258 
5 | 0093 | 0114 | 0119 0120 | 0120 | 0120 | 0120 } 0120 | 0120 
6 | 0044 | 0054 | 0056 | 0056 | 0056 | 0057 | 0057 | 0057 
7 | | 0021 | 0025 | 0026 | 0027 | 0027 | 0027 | 0027 
8 | | | 0010 | 0012 | 0012 | 0013 | 0013 | 0013 
9 | | | | 0005 | 0006 | 0006 | 0006 | 0006 
10 | 0002 | 0003 | 0003 | 0003 
11 | | | | | | 0001 | 0001 | 0001 
12 | 0000 | 0001 
13 | | 0000 
TABLE 9.—Values of 7;,(2) 
k |n=3 | n=4|n=5 | n=6 | n=7 | n=8 | n=9 |n=10/n=11|n=12 n=13|n=14 nz15 
1 | 1056 | 1250 | 1292 | 1301 | 1303 | 1304 | 1304 | 1304 | 1304 | 1304 | 1304 | 1304 | 1304 
2 | 3230 | 3750 | 3859 | 3883 | 3888 | 3890 | 3890 | 3890 | 3890 | 3890 | 3890 | 3890 | 3890 
3 1250 | 1493 | 1546 | 1558 | 1561 | 1561 | 1562 | 1562 | 1562 | 1562 | 1562 | 1562 
4 | 0539 | 0653 | 0678 | 0684 | 0685 | 0686 | 0686 | 0686 | 0686 | 0686 | 0686 
5 | 0245 | 0299 | 0311 | 0314 | 0314 oats 0314 0314 | 0314 | 0314 
6 | | 0114 | 0140 | 0145 | 0146 | 0147 | 0147 | 0147 | 0147 | 0147 
7 | | 0054 | 0066 0068 | 0069 | 0069 | 0069 | 0069 | 0069 
8 | | 0025 | 0031 | 0032 | 0033 | 0033 | 0033 | 0033 
9 0012 | 0015 | 0015 | 0015 | 0015 | 0015 
10 0006 | 0007 | 0007 | 0007 | 0007 
11 | | | 0003 | 0003 | 0003 | 0003 
12 0001 | 0002 | 0002 
| 
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—NOTES— 


A METHOD FOR THE SOLUTION OF CERTAIN NON-LINEAR 
PROBLEMS IN LEAST SQUARES* 


By KENNETH LEVENBERG! (Frankford Arsenal) 


The standard method for solving least squares problems which lead to non-linear 
normal equations depends upon a reduction of the residuals to linear form by first 
order Taylor approximations taken about an initial or trial solution for the parame- 
ters.? If the usual least squares procedure, performed with these linear approxima- 
tions, yields new values for the parameters which are not sufficiently close to the ini- 
tial values, the neglect of second and higher order terms may invalidate the process, 
and may actually give rise to a larger value of the sum of the squares of the residuals 
than that corresponding to the initial solution. This failure of the standard method 
to improve the initial solution has received some notice in statistical applications of 
least squares* and has been encountered rather frequently in connection with certain 
engineering applications involving the approximate representation of one function 
by another. The purpose of this article is to show how the problem may be solved 
by an extension of the standard method which insures improvement of the initial 
solution.* The process can also be used for solving non-linear simultaneous equations, 
in which case it may be considered an extension of Newton's method. 

Let the function to be approximated be h(x, y, z, - - - ), and let the approximating 
function be H(x, y, z,---;a,B8,y,+++),wherea,B,y, --- are the unknown param- 
eters. Then the residuals at the points, (x:, yi, 2i,°-+),7=1,2,°++-+,m, are 


fila, B, 4 eee ) = A(x, Vis Sin °° ° 5 ey B, wero ) - h( xi, Vi» = * ); (1) 


and the least squares criterion requires the minimization of 


s(a,B,y, °°: ) = Df. (2) 


(It is assumed that the weights of the residuals are unity. If not, consider the func- 


* This paper was read before the Annual Meeting of the American Mathematical Society in Chicago, 
Ill., on Nov. 26, 1943. Manuscript received Feb. 2, 1944. 

1 The writer wishes to thank Dr. J. G. Tappert, under whose direction the method of damped least 
squares was developed, and Dr. H. B. Curry, for valuable suggestions and guidance. 

2E, T. Whittaker and G. Robinson, The calculus of observations, Blackie and Son, London, 1937, 


p. 214. 
3 E. B. Wilson and R. R. Puffer, Least squares and laws of population growth, Proc. Amer. Acad. Arts 


and Sci. (Boston), 68, 285-382 (1933). 
4 Another extension of the standard method, which requires the use of second partial derivatives, 


is given by Wilson and Puffer (I.c.). 

A different kind of approach, not based upon the standard method, is given by Cauchy, Méthode 
générale pour la résolution des systémes d’ équations simultanées, C. R. Acad. Sci. Paris, 25, 536-538 (1847). 
See also a paper by H. B. Curry, not yet published, (abstract in Bull. Amer. Math. Soc., 49, 859 (1943), 


abstract No. 278). 
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tion f; to be the product of the residual and the square root of the corresponding 


weight.) Choosing an initial solution, po = (ao, Bo, Yo, - * - ), at which it is assumed 
that s does not have a stationary value, the first order Taylor expansions of the 
residuals are taken about fo, giving a set of linear approximations to the residuals, 


2 a Of; Of; Of; 
Jila, B, > sees ) = F(a, B, Y: aiid ) = fi(po) + —— Aa + oo AB + —— Ay + a (3) 
da ap ay 


where Aa=a—ao, AB=B—Be, +--+, and the partial derivatives are evaluated at po. 
Now, the standard method consists of minimizing 


S(a,.B,y°:) = DF (4) 
1 


by setting the partial derivatives of S with respect to the various parameters equal 


to zero, yielding the usual linear normal equations, 


1 0S 

ss -= [aa |Aa + [a6 |AB + lay |Ay +---+ [a0] = 0, 

1 0S (5) 
> ag 7 ala + [06]06 + [6vlav + --- + [60] = 0, 


where the notation [ ] is a symbol of summation, so that, e.g., 


l= Ef ol E(B), at E (CE), ate 


1 Ja 1 1 





However, as pointed out above, the values of the increments, Aa, AB, Ay, - - - , ob- 
tained by solving equations (5), may be so large in absolute value as to invalidate the 
approximations (3) so that the decrease in S may not correspond to a decrease in s. 

In such cases, it would seem advisable to limit or “damp” the absolute values of 
the increments of the parameters in order to improve the first order Taylor approxi- 
mations (3) and to minimize simultaneously the sum of the squares of the approximat- 
ing residuals (4) under these damped conditions. In order to make both the incre- 
ments and the residuals small in absolute value, the least squares idea can be em- 
ployed. The sum of the squares of both the residuals and the increments may be 
minimized. More precisely, the expression to be minimized will be 


S(a, B, y,°*+) = wS(a, B, y,-++*) + a(Aa)? + b(AB)? + c(Ay)? +--+: , (6) 


where a, 3b, ¢, are a system of positive constants or weighting factors expressing 
the relative importance of damping the different increments, and w is a positive 
quantity expressing the relative importance of the residuals and increments in this 
minimizing process. If we denote the point at which S takes its minimum, for any 
positive value of w, by pu=(aw, Bu, Yw, * * * ), and set 


O(a, B, y,-- +) = a(Aa)? + (AB)? + c(Ay)? +---, (7) 


it is seen, under the assumption that s is not stationary at fo, that 
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WS (Pw) < wWS(Pw) + O(pw) = S(pu) < S(po) = wS(po) + Q(po) = wS(po), 
whence S(pu) < S(po). (8) 


Also, denoting the standard least squares solution by p,, (the reason for the notation 
is discussed later), we have 


WS (Pw) + O(pu) = S(pw) < S(pu) = wS(px) + O(Px) < WS(pw) + O(px), 
whence O(pw) < Opa). (9) 


Inequality (8) shows that the minimization of (6) will diminish the sum of the squares 
of the approximating residuals, S, and (9) shows that the increments given by the 
standard least squares solution will be improved in the sense that the weighted sum 
of their squares, Q, will be reduced. That the sum of the squares of the true residuals, 
s, can be diminished, will be proved shortly. 

To minimize (6) and obtain p,, the partial derivatives of S with respect to the 
various parameters are put equal to zero, and we get 
os 0s as os 

-=4% --+ 2aAa = 0, — = w— + 20AB = 0, 

da Oa Op 0B 

When we divide through by 2w, and substitute the expressions for the partial deriva- 


tives of S from (5), the “damped normal equations” become 
[aa] + aw) Aa + [aB|AB + [ay]Ay + --- + [a0] = 0, 
[Ba|Aa + ([88] + bw-)AB + [By]Ay + --- + [80] = 0, (10) 


These equations are seen to be the same as the ordinary normal equations (5), except 
for the coefficients of the principal diagonal, which are increased by quantities pro- 
portional to the weighting factors a, b, c,- ++, respectively. Since the symmetry 
of the matrix of the coefficients of equations (5) is preserved, simplified methods of 
solution of linear simultaneous equations, which take full advantage of such sym- 
metry,> may be used to solve equations (10). It is to be noted that the standard 
method of least squares corresponds to w— ~, and is thus a special case of the method 
here given, which may be termed the method of “damped least squares.” 

If we denote the number of parameters by k, it is seen from the determinantal 


solution of equations (10) that, in the neighborhood of w=0, 





= [a0 w'-'hed +--+ 4+ --- 
Aa = ay, — ajo=- ——_— oe = [a0 Ja—' w t+-:: 
w *abe eee + > em 
day 
whence 7; - =— [a0 Ja - (11) 
aw w=0 


and similarly for the other parameters. Now 
ds( Pw) ds da ds dp 


dw da dw OB dw 


5 P. S. Dwyer, The solution of simultaneous equations, Psychometrika, 6, 101-129 (1941). 


1944] KENNETH LEVENBERG 167 
and, from the definition of the summation symbols, we find that the partial deriva- 
tives of s at po are given by 


0a 


Os Os 

— = 2]a0, — =? ae 13 
) el (13) 
Hence the substitution of (11) and (13) in (12) yields 


ds 
(=) = — 2{ [e0]}*a— + [po]}*o-1 + -- - }. (14) 
IW / w=0 
This derivative is negative since the partial derivatives in (13) are not all zero, by 
the assumption that s does not have a stationary value at po. Therefore, s(pw) is de- 
creasing at w=0, thus insuring that values of w can be found for which the sum of the 
squares of the true residuals (2) will be reduced. 

The best value of w to use may theoretically be determined directly by solving 


ds( pw) 
— = 0; (15) 
dw 
however, this equation is generally complex in practice. By writing 
; ds 
S(pw) = s(po) + w{ — ‘ (16) 
dw w=0 


and setting the left side of (16) equal to zero on the assumption that pp» was chosen 
so that the decreased value s(p.) will be small, the approximate formula, 


S(po) >5(Po) 


, ” seins ae (17) 
ds/dw .0  [a0}?a— + [60] + ie 


is obtained.® If necessary, this value may be improved by calculating s(p~) for several 
different trial values of w, so that an approximate minimum may be located graphi- 
cally. Experience with the method, especially in connection with fitting a particular 
function H(x, y, 3, ---;a,8,7,-- +), enables one to get an idea of the general order 
of magnitude of the best value of w so that very few trial values of w should suffice. 
If so desired, the improved set of values of the parameters may be further improved 
(if the true minimum has not already been reached), by a repetition of the process, 
considering this improved set as a new initial solution. 

So far, the weighting system a, b,c, - - - has been left arbitrary, the only restric- 
tion being that the weighting factors be positive. If we set the criterion that these 
factors be chosen so that the directional derivative of s, taken at w=0 along the curve 
=A», B=By, - ++, should have its minimum value, namely, the negative gradient, 


we have 


ds f da\? dB \? —1/2 as\2 Os\2 1/2 
a) on ee ee a 
dw dw dw Oa 0p 


where the derivatives are taken at w=0. Substitution of (14), (11), (13) in (18) 


gives us 


6 This type of approximation was used by Cauchy (l.c.). 
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{ [0 ]*a- + [g0}%-!+.--- } { a[0]2a-2 + [80]b-2+--- }—1/2 
= {[a0]*+ [60]? +--+}, (19) 


and this is satisfied when the factors a, 6, c,- + - are all equal. Without loss of gen- 
erality, they may be taken equal to unity. For this weighting system, the formation 
of the damped normal equations (10) may be thought of as being accomplished simply 
by the addition of a positive constant, 1/w, to the coefficients of the principal diagonal 
of the standard normal equations (5). Another weighting system which has been used 
successfully is, a= [aa], b=[88], - - - ; in this case the damped normal equations 
are formed by multiplying the principal diagonal coefficients of the standard normal 
equations by a constant greater than unity, 1+1/w. 

The nature of the damping which we have imposed upon the parameter variables 
can be given a simple geometric interpretation. For instance, if the unity weighting 
system is considered, the “overshooting” of the solution is prevented by damping the 
distance (& dimensional) from the initial solution point, since Q is then the square 
of this distance. By this restriction of k dimensional distance (which would appear 
to be a natural way to prevent overshooting), we are not obliged to decide on an ar- 
bitrary preassigned procedure restricting the variables individually, as is done, for 
example, by the method of Cauchy (l.c.). The greater freedom given the individual 
variables by the method of damped least squares may account for the fact that it 
has solved, with a comparatively rapid rate of convergence, types of problems which 
are of much greater complexity than those to which the principle of least squares is 


ordinarily applied. 


ON THE DEFLECTION OF A CANTILEVER BEAM* 
By H. J. BARTEN (Washington Navy Yard) 


In spring theory it is sometimes necessary to compute the deflection of a cantilever 
beam for which the squares of the first derivatives cannot be neglected as is done in 
classical beam theory. This problem is thus placed in the same category as the prob- 
lem of the elastica. 

The solution given in this note can be applied to a cantilever of any stiffness. The 
difference between the deflection as found by the classical beam theory and that 
found by the present method is, however, noticeable only in the case of beams of 
low stiffness. 

The clamped end of the beam is taken as the origin of coordinates and downward 
deflections are considered as positive. A point on the beam may be identified by four 
quantities of which only one is independent. These four quantities are the two rec- 
tangular coordinates x and y, the arc length s measured from the origin of coordinates, 
and the deflection angle 6 which is the angle between the tangent to the curve at the 
point under discussion and the horizontal. We may thus identify this point by the 
symbol (x, y, s, 8). The subscript Z is used to identify the value of these quantities 
at the free end of the beam. Before deflection a vertical load P is applied at the point 
(L, 0, L, 0). The beam has a uniform cross section of moment of inertia J and is com- 





* Received Feb. 21, 1944. 
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posed of a material whose modulus of elasticity is E. The problem is to find the de- 
flection of the end-point of the beam due to the vertical load P. 
The bending moment induced at the point (x, y, s, 6) by the vertical load P is 


M = P(x, — x). 
Therefore 
do/ds = a(x, — x), (1) 
where a=P/EI. Using the relation 
dé d@ dx dé 
— = — — = cos# —» 
ds dx ds dx 
we obtain 
f 0s 6d@ = ff aes — x)dx 
or 
sin 0 = a(xypx — $2x°) +C. (2) 


The boundary condition at the clamped end of the beam, namely, @=0 when x =0, 


reduces Eq. (2) to 
sin 06 = a(xypx — 3x?). (3) 


Thus . 
sin 6, = $axz. (4) 


Combining the latter expression and Eq. (3) we obtain 
sin 0, — sin @ = a(x, — x). (5) 


Thus 
ez — = [20+ (sin 0, — sin)". 


Substituting this expression in.» Eq. (1), we obtain 








dé dé dy } dé : ; 
— = — — = sin @ — = [2a (sin 6, — sin 6) |'/2, 
ds dyds dy 
or 
° sin 6 d@ 
7 ones rage te 
0 [2a(sin 0; — sin @) }'/2 
There fore 
oL sin 6 d@ 
VL -{ [7 (3 : “i (6) 
0 [2a(sin 6, — sin 6) ]'/? 


With the transformation 


(5-5) -o(G-F)see= ts 
cos§ — — —] = cost — — — one = 28m @, 
“a 3 4 2 ' 


uf r/2 (2k? sin? ¢ —1)d@ 7) 
yp = oll? —y 7 
™ 3 (1 am Rk? sin? o)'!? 


Eq. (6) becomes 





where 
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. cos 7/4 T Or 
sin 6 = ——— —9 k = cos {[ — —- ‘ 
4 ? 


Eq. (7) is a combination of incomplete and complete elliptic integrals! and may be 
written 


VL = a~/2 F(R) — F(k, 6) — 2E(k) + 2E(R, 6) |, (8) 


where F(R) and E(k) are the first and second complete elliptic integrals respectively 
and F(k, 6) and E(k, 6) are the first and second incomplete elliptic integrals respec- 


tively. 
As Eq. (8) stands it is useless unless we find 6, as a function of a and L. This rela- 
tionship may be obtained in the following manner. From Eq. (1) we get 


L 
Or = f Gin = x)ds. 
0 


Integrating by parts we obtain 


ry L dx L 
0; = asdz = as—ds = as cos 6 ds. 
0 0 ds 0 


Differentiating this latter integral with respect to its upper limit, we have 
d6,/dL = aL cos @,. 


The solution to this differential equation is 


; aL? 
sin 6, = tanh —* 


(9) 


This completes the solution to the problem. 

In order to compare our results with those of Gross and Lehr? we must express 
our solution in the same dimensionless factors that they employed. By dividing the 
actual deflection of the beam by the “small deflection” aZ*/3 they obtain a deflection 
factor which is a function of the dimensionless quantity aL. We shall call this deflec- 


tion factor F,. Thus, from Eq. (8) 
3yz 
i aL? 
In order to find the maximum bending stress at the clamped end of the beam we 
must know the length of the moment arm x,. Combining Eqs. (4) and (9) we find that 


F = 3(aL*)-*/2[F(k) — F(k, 6) — 2E(k) + 2E(k, 6) |. (10) 


2 aL? ° 
x; = — tanh —- (11) 
a Z 


Gross and Lehr use the dimensionless contraction factor x,/Z an an aid in find- 
ing xz. We shall define this factor as F,. Thus 


1 Jahnke and Emde, Funktionentafeln mit Formeln und Kurven, Dover Publications, 1943. 
2 Gross and Lehr, Die Federn, V. D. I. Verlag, 1938. 
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2 aL? 
F? = — tanh —-- (12) 
aL? 2 


Computations show that Gross and Lehr’s values of F, have a constantly increas- 
ing error which deviates about 4% from our results when aL?=1. 
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The two factors F, and F, are very important to the designer. For this reason 
curves of these two factors with aL? as the independent variable are given in Fig. 1. 
The values of F, were computed from Jahnke and Emde. 


ON WAVES IN BENT PIPES* 
By S. A. SCHELKUNOFF (Bell Telephone Laboratories) 


In a recent issue of this QUARTERLY,! Karlem Riess obtained expressions for the 
fields of electromagnetic waves in bent pipes of rectangular cross section by the 
perturbation method. While it is true that in a bent pipe the waves cannot be classi- 
fied into transverse electric and transverse magnetic types because in general both 
E and H have components in the direction of wave propagation, a different classifica- 
tion into two types is possible. This permits another method which yields the general 
solution in terms of Bessel functions. 

In the one wave type, the plane of the electric ellipse is normal to the axis of 
bending (the Y-axis in Figure 1, p. 329 of Riess’ paper); these waves have been called 
electrically oriented (EOm,» Wave type) and the fields of these waves are obtainable from 
H, which may be expressed as the product of Bessel and sine (or cosine) functions. 

* Received Feb. 18, 1944. 

1 Vol. 1, No. 4, pp. 328-333. 
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In the other wave type, the plane of the magnetic ellipse is normal to the axis of 
bending; these waves are magnetically oriented (MOm,, wave type) and their fields are 
obtainable from E,. In each case the order of Bessel functions is equal to the angular 


phase constant. 
For a bent pipe formed by the intersection of two concentric spheres and two 


coaxial cones emerging from the center there is also a solution in terms of known func- 
tions. In one wave type, EO, ,, type, the plane of the electric ellipse is normal to the 
radius; in the other, MOn,, type, the plane of the magnetic ellipse is normal to the 
radius. The fields of EO-waves are calculable from H, and the fields of 1/O-waves 
from E,; H, and E, themselves can be expressed in terms of Bessel and Legendre 
functions. These waves may be called spherically oriented in order to distinguish them 
from the plane oriented waves described earlier. The letters S and P in front of EO 
and MO may be conveniently used in the abbreviations. 


CORRECTIONS TO MY PAPER 


A STRAIN ENERGY DERIVATION OF THE TORSIONAL- 
FLEXURAL BUCKLING LOADS OF STRAIGHT COLUMNS 
OF THIN-WALLED OPEN SECTIONS 


QUARTERLY OF APPLIED MATHEMATICS, 1, 341-345 (1944). 


- 


By 
N. J. HOFF 
In the last term of the right hand side member of Eq. (3) on page 343, should 


be raised to the second power and not to the fourth power. 
The following equation defining T should be added: 


T = (1/p?)(m?R + GC). 
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